Modelling Wave Propagation in Two-dimensional Structures
using a Wave/Finite Element Technique
E. Manconi and B.R. Mace
ISVR Technical Memorandum 966
February 2007
SCIENTIFIC PUBLICATIONS BY THE ISVR
Technical Reports are published to promote timely dissemination of research results
by ISVR personnel. This medium permits more detailed presentation than is usually
acceptable for scientific journals. Responsibility for both the content and any
opinions expressed rests entirely with the author(s).
Technical Memoranda are produced to enable the early or preliminary release of
information by ISVR personnel where such release is deemed to the appropriate.
Information contained in these memoranda may be incomplete, or form part of a
continuing programme; this should be borne in mind when using or quoting from
these documents.
Contract Reports are produced to record the results of scientific work carried out for
sponsors, under contract. The ISVR treats these reports as confidential to sponsors
and does not make them available for general circulation. Individual sponsors may,
however, authorize subsequent release of the material.
COPYRIGHT NOTICE
(c) ISVR University of Southampton
All rights reserved.
ISVR authorises you to view and download the Materials at this Web site ("Site")
only for your personal, non-commercial use. This authorization is not a transfer of
title in the Materials and copies of the Materials and is subject to the following
restrictions: 1) you must retain, on all copies of the Materials downloaded, all
copyright and other proprietary notices contained in the Materials; 2) you may not
modify the Materials in any way or reproduce or publicly display, perform, or
distribute or otherwise use them for any public or commercial purpose; and 3) you
must not transfer the Materials to any other person unless you give them notice of,
and they agree to accept, the obligations arising under these terms and conditions of
use. You agree to abide by all additional restrictions displayed on the Site as it may
be updated from time to time. This Site, including all Materials, is protected by
worldwide copyright laws and treaty provisions. You agree to comply with all
copyright laws worldwide in your use of this Site and to prevent any unauthorised
copying of the Materials.
UNIVERSITY OF SOUTHAMPTON
INSTITUTE OF SOUND AND VIBRATION RESEARCH
DYNAMICS GROUP
Modelling Wave Propagation in Two-dimensional
Structures Using a Wave/Finite Element Technique
by
E. Manconi and B.R. Mace
ISVR Technical Memorandum No: 966
February 2007
Authorised for issue by
Professor M.J. Brennan
Group Chairman
© Institute of Sound & Vibration Research
Abstract
The purpose of this work is to present a general method for the numerical analysis of wave
propagation in 2-dimensional structures by the use of a finite element method (FEM). The
method involves typically just one finite element to which periodicity conditions are applied
instead of modelling the whole structure, thus reducing drastically the cost of calculations. The
mass and stiffness matrices are found using conventional FE software. The low order dynamic
stiffness matrix so obtained is post-processed and the wavenumbers and the frequencies then
follow from various resulting eigenproblems.
The method is described and numerical examples given. These include isotropic and orthotropic
plates, isotropic cylindrical shells and the more complex case of sandwich cylindrical shells for
which analytical solutions are not available. The last two cases are studied by postprocessing
an ANSYS FE model.
The main advantage of the technique is its flexibility since standard FE routines can be used
and therefore a wide range of structural configurations can be easily analysed. Moreover the
propagation constants for plane harmonic waves can be easily predicted for different propagation
directions along the structure. The method is seen to give accurate predictions at negligible
computational cost.
1
Contents
1 Introduction
3
2 The
2.1
2.2
2.3
2.4
2.5
Wave–FE method
Finite element analysis of 2-dimensional structures . . .
Wave propagation in 2-dimensional uniform structures
Periodicity conditions . . . . . . . . . . . . . . . . . . .
Real propagation constants . . . . . . . . . . . . . . . .
Complex propagation constants . . . . . . . . . . . . .
3 Plane wave propagation in plates
3.1 Isotropic plate . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Real–valued dispersion curves and wave forms
3.1.2 Accuracy . . . . . . . . . . . . . . . . . . . . .
3.1.3 Complex–valued propagation constants . . . .
3.1.4 Sensitivity analysis of eigenvalues . . . . . . .
3.1.5 Sensitivity analysis of the wavenumbers . . . .
3.2 Orthotropic plate . . . . . . . . . . . . . . . . . . . .
3.2.1 Real–valued dispersion curves . . . . . . . . .
3.2.2 Complex–valued propagation constants . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
6
6
7
9
9
9
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
12
16
18
19
20
23
23
26
4 Wave propagation in curved and cylindrical shells
27
4.1 Isotropic shell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Sandwich shell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5 Concluding remarks
40
2
1. Introduction
One of the most general approximate methods for solving problem concerning the dynamics
of continuous structures is the finite element method (FEM). In its “classical” modal/dynamic
application, this method is used to obtain information about the vibrational behaviour from
the model of the whole structure in the low frequency range. However, in many engineering
applications, high frequency vibrations become significant, in particular where sound transmission have to be considered. At high frequencies, when the wavelengths become very small, the
use of this approach may become inappropriate because the size of the FE model that can
predict response accurately becomes too large and therefore computationally very expensive.
Furthermore, numerical inaccuracies arise at higher frequencies. One of the main reasons is the
dispersion error involved in the FE approximation. The dispersion error refers to the difference
between the numerical estimates of the wavenumbers and the true values. Since this error
increases with the frequency, in order to keep it within “acceptable” limits, the FE model size
has to be increased leading to a drastic increase in the computer time and storage. Therefore,
at high frequencies, a detailed FE model may be very difficult to solve or perhaps impossible
sometimes. One alternative technique to FEM is Statistical Energy Analysis (SEA). It is mainly
used for structures that can be considered as an assemblage of subsystems. The application of
SEA requires the knowledge of the modal density and of the relation between modal energy and
vibration velocity for each subsystem. However, while for simple structures these characteristics are well known, for complex structures the analysis is often very difficult. Computational
techniques such as spectral FEM or the technique proposed in the present report can be used
to evaluate these quantities.
Another approach to obtain information about vibration of structures is the Dynamic Stiffness
Method, DSM, sometimes also termed Spectral Element Method when it is applied to waveguides, [6, 7]. The method divides the structure into simple elements whose degrees of freedom
are defined at certain points called nodes. The key of this technique is the establishment of a
dynamic stiffness matrix in the frequency domain to relate nodal responses x and the forces F,
i.e.
D(ω)x = F.
Only one element has to be analysed and a global dynamic stiffness matrix can be subsequently
defined by assembly in a similar manner to FEM. Since the matrix D is obtained from the exact
solutions of the differential equation of motion, the method describes exactly the behaviour of
the element but it can only be applied to relatively simple cases. Flotow, [8], and Beale and
Accorsi, [9], applied these ideas to study the vibrations of network structures composed of simple
beam elements. They analysed the dynamics of the structure in terms of wave propagation in
a single beam element. Hence, they solved for the dynamics of the whole structures taking
into account the boundary conditions between every single element and those for the whole
structure. Although the method that they proposed reduces the calculation cost greatly, it is
developed for a particular class of structures like frame structures.
A large amount of research has focused on the question of extending the range of frequency
application of the FEM. In this context one method for structural dynamics and acoustic
applications is the spectral finite element method (SFEM). The SFEM strategy applies to
structures that present symmetry in one direction, called a waveguide. The name is due to
the fact that they allow waves to propagate in the direction of the axis of symmetry of the
structure. As it will be described later they are a particular case of periodic structures.
3
The basic steps of the SFEM are described in [1]. The SFEM uses the FEM to describe the
cross–section motion. At the same time the nodal displacements along the axis with which
the structure is aligned are in effect assumed to be harmonic functions of space. For harmonic
vibrations, the solutions to the wave equation are given by a polynomial eigenvalue problem.
These basic assumptions have been used by several authors to describe wave propagation in
different kinds of waveguide structures. Some examples include the paper by Gavric who
applied the method to a railway track, [2], Finnveden, who studied rib-stiffened structures and
fluid-filled pipes, [3, 4] and Shorter who developed SFE for viscoelastic laminates, [5].
Other authors have proposed a different technique to analyse wave propagation in periodic
structures. The term “periodic structures” indicates structures which exhibit characteristics
that repeat periodically in one, two or three-dimensions. In other words, a periodic structure
can be considered as an assemblage of identical elements, called cells or periods, which are
coupled to each other on all sides and corners by identical junctions. This characteristic is
indeed observable in many real systems. Examples include railway tracks, flat or curved panels
regularly supported, such as stringer stiffened panels, fluid filled pipes with flanges, acoustical
ducts, rail structures, car tyres, composite plates or shells etc. The use of the structures periodicity can be exploited to solve dynamic problem with reduced computational cost. Many
previous works have concerned periodic structures. One of the classical books where the mathematics of wave propagation in periodic structures has been discussed is that of Brillouin [10].
More recent works, [11–15], have studied the problem of defining the dynamical behaviour of
periodic structures considering the properties of its characteristic waves. These works consider
structures with cross–sectional properties that are uniform along one direction. Much attention
is focused on uniform plates and shells with identical stiffeners at regular intervals. In particular in [12] a finite element approach is used for obtaining the dispersion curve of a periodic
structure. The main idea is to relate the displacements, q, and the forces, f, at the left and
the right hand–side of a periodic element of the structure by the propagation constant λ in the
following way:
qR = λqL , fR = −λfL ,
√
where λ = eiµ , µ is the propagation constant and i = −1.
Following this idea, in [16–18], an estimate of the wave motion along a uniform waveguide
is carried out by considering a transfer matrix to relate the displacements and the forces on
both left and right hand–sides of a section of the waveguide. The transfer matrix is the result
of a manipulation of the mass and stiffness matrices obtained from a conventional FE model
of a small section of the waveguide. In this Wave–FE (WFE) approach, wave propagation is
described in terms of eigenvalues and eigenvectors of this matrix and the dispersion relation
is obtained via a standard eigenvalue problem. In particular, in [16], after reformulating the
dynamic stiffness matrix of one cell of a waveguide in terms of wave propagation, the dynamic
stiffness matrix of the whole structure is built using periodic structure theory. The response of
the whole structure can then be evaluated simply. The two examples provided in [16], beams
and plates excited by a point force, have shown that the accuracy of the method is good when
common requirements for FE discretisation are satisfied, that is, the length of the element is
small enough compared to the wavelength. In the same way as the SFEM, the advantage of
this strategy is that only one section of the structure has to be meshed and solved, reducing
drastically the cost of calculation. However, unlike the SFEM, the method described earlier
exploits the high geometrical and material flexibility of standard commercial FE packages to
extract the stiffness, mass and damping matrices. This is a great advantage since no new
elements have to be developed for each application. Furthermore, the SFEM method cannot
4
be applied to 2-dimensional structure.
The purpose of this work is to extend the WFE approach to 2-dimensional structures. The
advantages in extending the method to the 2-dimensional case consists mostly of the fact that
only one finite element can be used and that the propagation constants of plane harmonic
waves can be easily evaluated for different propagation direction across the structure. An
important reference for this work is [19], where FEM is applied to model generic periodic
structures. Herein, a general method is presented that enables vibrational waves propagating
in a periodic 2-dimensional structure to be analysed by the use of a finite element method.
Similar to the approach presented in [16–18], this method relates the right-left and the topbottom displacements and generalised forces of one rectangular element by
qR = λx qL ,
fR = −λx fL ,
qT = λy qB ,
fT = −λy fB ,
where λx = eiµx , λy = eiµy and µx , µy are the propagation constants of a plane harmonic wave in
a 2-dimensional geometry. Making use of the FEM, the mass and stiffness matrices are obtained
and they are subsequently post-processed in order to obtain a problem whose eigenvalues and
eigenvectors describe the wave propagation characteristics. Special attention is given to this
eigenvalue problem. In fact if µx and µy are assigned, the frequencies of wave propagation
can be obtained from a standard eigenvalue problem. On the other hand if the propagation
frequency is prescribed it yields a polynomial or a transcendental eigenvalue problem. The
real dispersion curves can be easily obtained by solving the first case. The solutions to the
transcendental eigenproblem are less obvious nevertheless, for the more general case of complex
propagation constants, they represent decaying wave motion at a given frequency. For the
first kind of problem efficient algorithms, which include error and condition estimates, are well
known while the last problem, which is more general, is less obvious and different methods to
solve it are analysed.
This memorandum is organised in three parts. The first part concerns the description of the
general method by which the wavenumbers for a 2-dimensional structure can be predicted from
a FE model. The second and third parts deal with the application of the method to several
examples. First an isotropic plate is considered as an illustrative example. Another example
concerning an orthotropic plate made of a Glass/Epoxy composite material is discussed. In
the third part cylindrical shells structures are discussed. In particular an example of a curved
sandwich panel, for which analytical solutions are not available, is provided. The dispersion
relations, the complex propagation constants and wave-forms are given. A sensitivity analysis
is also performed in order to identify those results which truly represent the behaviour of the
continuous structure. Finally the results are discussed and conclusions are drawn.
5
2. The Wave–FE method
In this section a general theoretical approach to apply the wave-FE method to a 2-dimensional
periodic structure is described. The main idea is to use a wave approach instead of a traditional modal one to describe the vibration behaviour of a uniform structure modelled by a FE
technique.
2.1. Finite element analysis of 2-dimensional structures
A two dimensional periodic structure can be considered as an assembly of identical elements
which are coupled to its neighbours on all sides and corners. Every element, or cell, of the
system can be indicated by two numbers (ni , nj ) which defines its position along the directions
d1 and d2 parallel to the direction of the system periodicity as shown in Figure 1. The dynamical
behaviour of such structures can be studied in a simple way if the behaviour of only one cell is
analysed. The general procedure for FE discretisation of continuous structures is well known
and therefore will not be described here. The cell can be in general represented by a model
d2
n1 , n2 + 1
n1 − 1, n2
n1 , n2
n1 + 1, n2
n1 , n2 − 1
d1
Figure 1: Uniform 2-dimensional structure.
with internal and boundary degrees of freedom q. The only condition required for the following
analysis is that the cell is meshed with an equal number of nodes at its left–right and top–
bottom sides. Another requirement for the mesh is a small spatial discretisation compared
with the wavelength of the problem. As a rule of thumb, 6-10 elements per wavelength have to
be used to obtain accurate results.
Using the stiffness and mass matrices obtained from the FE model, the equation of motion for
the cell for time harmonic motion is
K − ω 2 M q = F,
(1)
where ω is the angular frequency, K and M are the stiffness and mass matrices, F is the loading
vector containing the nodal loads and q is the vector containing the nodal displacements. It is
worth noting that the loading vector F is responsible for transmitting the wave motion from
one element to the next and therefore, even if free wave motion is considered and so no external
loads exist, it is not zero.
6
z
y
3
4
Ly
1
x
2
Lx
Figure 2: Rectangular 2-dimensional element.
Now consider a uniform structure as a special case of a periodic structure. The structure
is meshed using conventional FEs and rectangular elements. The smallest cell that can be
considered in the analysis of the uniform structure is just one FE. It will be shown that, by
appropriately formulating the problem, this element is sufficient to estimate the dispersion
relations for the whole system. Figure 2 shows a 2-dimensional rectangular finite element. The
Cartesian axes x and y are chosen as the preferred directions d1 , d2 respectively. The element
is a rectangular four node element and x, y and z are the global coordinates. The nodes are
numbered as in Figure 2. The displacement vector for the element is q = [q1 , q2 , q3 , q4 ]T , where
qi contains the nodal degrees of freedom.
2.2. Wave propagation in 2-dimensional uniform structures
Plane waves propagating in the 2-dimensional structures are considered. Let φ be a generic
disturbance. It will propagate as a harmonic free wave if the physical problem admits a solution
of the type
φ = Aei(ωt−k·r) ,
(2)
where A is a constant, r is a vector with components x and y and k is a vector which represents
the wavenumbers for a wave travelling in the θ direction with respect to the x axis. In the
discrete model of the structure shown in Figure 1, at the points (ni Lx , nj Ly ),
φni ,nj = Aei(ωt−ni kx Lx −nj ky Ly )
(3)
where kx = k cos(θ), ky = k sin(θ), represent the components of the wavenumber k along the x
and y axes and Lx and Ly represent the element lengths. In the following analysis µx = kx Lx
and µy = ky Ly indicate the propagation constants for the structure while λx and λy equal e−iµx
and e−iµy respectively.
The propagation constants are usually complex numbers. The real parts of µx and µy represent
the change in phase in the wave motion from one cell to the next, while the imaginary parts of
µx and µy represent the wave attenuation along the x and y directions.
For the general case of a system discretised by a certain number of cells, Figure 1, the disturbance φ propagates through the structure such that
φn1 +1,n2 = λx φn1 ,n2 ;
φn1 ,n2 +1 = λy φn1 ,n2 ;
7
φn1 +1,n2 +1 = λx λy φn1 ,n2 .
(4)
Considering equation (4) for the simple case of one finite element as shown in Figure 2, we can
write the following relations for the nodal displacements:
q2 = λx q1 ;
q3 = λy q1 ;
q4 = λx λy q1 .
(5)
Equations (5) can be rearranged to give
q = Qr q1 ,
(6)
q = [q1 q2 q3 q4 ],
(7)
Qr = [I λx I λy I λx λy I]T .
(8)
where
and
Since equilibrium requires that the sum of the nodal forces at each node is zero, at node 1 then
F1
F2
−1
−1
(9)
I −λ−1
x I −λy I (λx λy ) I F = 0.
3
F4
Hence, premultiplying both side of equation (1) by
−1
−1
Ql = [I λ−1
x I λy I (λx λy ) I],
(10)
the equation of free wave motion takes the form:
[K(µx , µy ) − ω 2 M(µx , µy )]q1 = 0,
where
(11)
K = Ql KQr ,
(12)
M = Ql MQr
are the reduced stiffness and mass matrices. K and M are n×n matrices, where n is the number
of nodal degrees of freedom. The elements of these matrices are functions of the propagation
constants and of the elastic and inertial characteristics of the structure.
Equation (11) represents an eigenvalue problem in ω 2 for given values of µx and µy . This
equation can also be seen as an eigenvalue problem in µx and µy for a given ω and a given
propagation direction θ. This is a more general case but less easy to solve since it represents a
transcendental eigenvalue problem.
It can be seen from equation (11) that the mathematical formulation of the method is fairly
simple and that standard FE packages can be used to obtain the mass and stiffness matrices of
one period of the structures. This is a great advantage because there is no need to create FE
elements “ad hoc” every time that a new structure has to be studied.
8
2.3. Periodicity conditions
Equation (11) yields the same values of q and ω for (µx , µy ) or (µx + 2m1 π, µy + 2m2 π) where
m1 and m2 are integers, that is the propagation disturbance and the frequency are periodic
functions of the propagation constants. This is a consequence of considering a discretised
structure instead of a continuous one and represents a substantial difference with the standard
continuum approach. In fact, for the continuous model it is implicitly supposed that a certain
propagation characteristic could be “measured” between particles. Therefore in the continuous
model the observable wavelengths will be all the ones included in the interval 0 ≤ λ ≤ ∞.
On the other hand, since the FE is a discrete model of the structure, if the distance between
the nodes along x and y in one element are Lx and Ly , the model will allow any wavelength
components along x and y such that λx ∈ [0, 2Lx ], λy ∈ [0, 2Ly ], or, in term of the propagation
constants, µx ∈ [0, 2π], µy ∈ [0, 2π]. Since a wave propagates in the positive and negative
direction in the same way, Re(µx ) ∈ [−π, π], Re(µy ) ∈ [−π, π] are taken as the most convenient
interval in which to examine the variation of ω with respect to the variation of the propagation
constants. This restriction is not so arbitrary as it appears since it contains a complete period
of the frequency and avoid ambiguity in the wavenumbers at the same time. Other “intervals”
of the same kind can be found but are not of interest for the present work. Detailed discussion
about the periodicity conditions can be found in [10].
2.4. Real propagation constants
In this section we consider waves that propagate in the 2-dimensional structure without attenuation. This means that the propagation constants µx and µy are real quantities and for freely
propagating waves in an undamped structure
|λx | = 1,
|λy | = 1.
If the wave propagation constants are known, the problem of finding the propagating frequencies
mathematically means solving the standard eigenvalue problem for ω 2 , equation (11).
For real values of µx and µy , Qr = Q∗l , where ∗ indicates transpose conjugation. Due to the fact
that the stiffness and mass matrices in equation (1) are real, symmetric and positive definite
matrices, then for real propagation constants
∗
K = Qr KQl = Q∗l KQl = [Q∗l KQl ]∗ = K
(13)
∗
M = Qr MQl = Q∗l MQl = [Q∗l MQl ]∗ = M .
Hence, the dynamic matrices K and M are positive definite Hermitian matrices. For any given
values of the propagation constants there will therefore be n real positive eigenvalues for which
wave propagation is possible. This number is equal to the number of degrees of freedom per
node, n. It should be noted that, there is also a pair of possible waves travelling in the opposite
direction: if (µx , µy ) is a solution so is (−µx , −µy ).
2.5. Complex propagation constants
Unlike the previous analysis, equation (11) can be formulated also to give the propagation
constants for a given value of ω and θ. Two approaches are discussed to obtain the propagation
9
constants: with the first one the problem is reduced to a standard polynomial eigenvalue
problem while with the second the eigenproblem is a transcendental eigenvalue problem. For
the first kind of problem efficient algorithms, which include error and condition estimates, are
well known while the second problem, which is more general, is less common and different
methods can be used to solve it.
By partitioning the dynamic stiffness matrix of equation (1) as
D11 D12 D13 D14
D21 D22 D23 D24
D =
D31 D32 D33 D34 ,
D41 D42 D43 D44
equation (11) can be rewritten in the form
[(D11 + D22 + D33 + D44 ) + (D12 + D34 )λx + (D21 + D43 )λ−1
x +
−1 −1
+(D13 + D24 )λy + (D31 + D42 )λ−1
y + D14 λx λy + D41 λx λy +
(14)
−1
+D32 λx λ−1
y + D23 λx λy ]q = 0.
Considering the transpose and the transpose conjugate of equation (14) and exploiting the
symmetry of the dynamic stiffness matrix, it can be easily proved that the eigenvalues of this
−1 ∗
∗
−1
−1 ∗
problem occur in pairs (λx , λ∗x ), (λ−1
x , λx ) and (λy , λy ), (λy , λy ).
If the ratio µy /µx is a rational number, that is µy /µx = m2 /m1 , where m1 and m2 are integers
with no common divisors, we can reformulate the propagation constants as µx = m1 σ, µy =
m2 σ. Since λx = e−iµx and λy = e−iµy , equation (14) can be rearranged in this case as
[A0 + A1 γ m1 + A2 γ m2 + A3 γ 2m1 + A4 γ 2m2 + A5 γ m1 +m2 +
(15)
+A6 γ 2m1 +m2 + A7 γ m1 +2m2 + A8 γ 2m1 +2m2 ]q = 0,
where γ = e−iσ and the matrices Aj are of the same order as the partitions of D. Equation
(15) represents a standard polynomial eigenvalue problem of order n(2m1 + 2m2 ) and can be
solved by efficient algorithms. The eigenvalues of equation (15) are transcendental functions
in the complex variable σ. In this particular formulation, the problem admits an infinity of
periodic solutions since, if σ is a solution, σ ± 2pπ, p integers, is a solution as well. In terms of
the propagation constants µx and µy , the correspondent periods for the roots of equation (15)
will be 2m1 π and 2m2 π.
In order to consider a general way to solve equation (14) for any possible (irrational) value of
µy /µx we can write, for waves propagating at an angle θ with respect to the x axis
λx = e−ikLx cos(θ) ,
(16)
λy = e−ikLy sin(θ) .
This leads to a general formulation of the problem
A(e−ikLx cos(θ) , e−ikLy sin(θ) )q = 0,
10
(17)
where the elements of the matrix A are transcendental functions of the complex variable k for
a given value of θ.
The analytical functions involved in equation (17), sums and products of exponential functions,
are continuous and continuously differentiable with respect to the variable k. Therefore, the
most natural technique to solve (17) is to reduce the problem to the transcendental equation
g(k) = |A(k)| = 0 and solve it by Newton’s method. Newton’s method seeks the solutions of
the nonlinear equation by solving iteratively a sequence of linear problems of the form
g(ki−1 ) + ri g ′ (ki−1 ) = 0
(18)
ki = ki−1 + ri .
Sufficient conditions for the existence of the solution and the convergence of Newton’s method
are given by the Kantorovich Theorem, see [20].
A similar way to solve the problem is to use the Newton’s eigenvalue iteration method [21], [22].
This method extends Newton’s method to the matrix A(k) in the following way. Given
B=−
then solve for
dA(k)
,
dk
[A(ki−1 ) − ri (B(ki−1 )]qi = 0,
(19)
ki = ki−1 + min(ri ),
where min means the minimum absolute value of the eigenvalues ri . The main difference with
the same method applied to the one-dimensional case is that the values of ri are determined
by solving an eigenvalue problem and the approach is more general. Although both problems
consider first order approximations, higher order approximations can be included by considering
higher order terms in the series expansion of either the matrix A or of the function g.
An alternative choice seeks the complex roots of the equation |A| = 0 using a variant of the
Powell Method, [23]. This algorithm is implemented in the function fsolve in the Optimization
Toolbox of MATLAB as a default choice.
The solutions obtained by applying the three different methods are virtually identical but the
function fsolve is computationally more efficient than the other two. Despite the fact that the
solution space is more sensitive to Re(k) than Im(k), the above algorithms used have good
guaranteed convergence within the closed intervals used to find the solutions. Moreover, the
objective function |A| at the estimated solutions k was found to always be very “small”.
However two important aspects are still open to discussion: the number of solutions and the
interval in which all these solutions are defined. The choice of this interval is a more intersting/relevant problem as all the iterative root–finding methods depend critically on the initial
estimates.
Alternative numerical approaches for finding complex roots of the transcendental equation (17)
may be the Interval Newton method, [24], [25], the contour integration method, [26] or Muller’s
method, [27]. The evaluation and comparison of these methods is object of future work.
11
3. Plane wave propagation in plates
In this section flexural wave motion in a 2–dimensional thin plate is investigated. First an
uniform isotropic plate is analysed as an illustrative example. The second part concerns the
more complicated and practical case of an orthotropic plate made of Glass/Epoxy composite
material. The real–valued dispersion curves, the complex propagation constants and accuracy
and sensitivity analysis of the results are given. Although the complex propagation constants
for given frequency and for given propagation direction have been obtained, the procedure that
enables the complete complex–valued dispersion curves to be fully evaluated is in progress.
As stated before, only one single finite element is analysed, Figure 2. For a plate in bending
vibration, the element has three degrees of freedom at each node: translation in the z direction
and rotations about the x and y axes. The shape function assumed for this element is a complete
cubic to which the two quadratic terms x3 y and xy 3 have been added. For more detail see [28].
3.1. Isotropic plate
A steel plate of thickness h located in the (x, y) plane is considered. In the numerical examples
below the material properties are taken to be: Young’s modulus E = 19.2·1010 N/m2 , Poisson’s
ratio ν = 0.3, density ρ = 7800 kg/m3 .
The governing equation for the out–of–plane motion w of the thin plate is:
4
∂ w
∂2w
∂4w
∂4w
D
=
−ρh
+
2
+
,
(20)
∂x4
∂x2 ∂y 2
∂y 4
∂t2
where
Eh3
,
D=
12(1 − ν 2 )
is the flexural rigidity. A nondimensional frequency is defined as
p
Ω = ωL2x ρh/D.
If nothing different is specified the plate element is assumed to be a square element of length
L. Harmonic wave solutions of the form
w = Aeiωt−ikx x−iky y
(21)
are sought. Substituting (21) into equation (20), gives
D(kx2 + ky2 )2 = ρhω 2 .
(22)
Since kx = k cos(θ), ky = k sin(θ), and hence k = kx2 + ky2 , the analytical plate free wavenumber
is defined as
r
√ 4 ρh
.
(23)
k= w
D
3.1.1. Real–valued dispersion curves and wave forms
Figure 3 shows the plots of the nondimensional (Ω, µx , µy ) surfaces, where µx = kx Lx and
µy = ky Ly for real kx and ky . As expected the surfaces are symmetric with respect to the
12
18
16
14
12
Ω1
10
8
6
4
2
0
1
1
0.5
0.5
0
0
−0.5
−0.5
−1
µy/π
−1
µx/π
50
40
Ω2
30
20
10
0
1
1
0.5
0.5
0
0
−0.5
−0.5
−1
µy/π
−1
µx/π
60
55
50
45
Ω3
40
35
30
25
20
15
1
1
0.5
0.5
0
0
−0.5
µ /π
−0.5
−1
−1
y
µ /π
x
Figure 3: Isotropic plate. Propagation surfaces.
13
µx = 0 and µy = 0 planes. This symmetry is due to the isotropy of the plate and due to the
fact that a wave propagates in the same way in the positive or in the negative direction. The
wave motion is possible only at the frequencies specified by the surfaces in Figure 3 since at
other frequencies the wave decays exponentially. Hence we call these surfaces as passing bands
according to Brillouin, [10].
In order to highlight some important properties of the wave-FE model, Figure 4 shows the
numerical and analytical dispersion curves for a wave propagating along the x direction. It
Ω
60
Numerical. First branch
Numerical. Second branch
Numerical. Third branch
Analytical
Ωc
Third
propagating
band
50
40
Second
propagating
band
30
20
Ωb
Ωa
First
propagating
band
0
−3
Accurate
−2
−1
−0.62
0
µx/π
0.62
1
2
3
Figure 4: Real–valued dispersion curves for θ = 0. Numerical and analytical curves.
can be seen that the frequency is periodic with a period 2π. Therefore, as explained in the
section 2.3, it is sufficient to evaluate the variation of the frequency within the interval [−π, π].
As expected three passing bands are obtained for the plate because there are three DOFs per
node. However, only the first, which corresponds to the results obtained from the analytical
formulation, exists for the continuous plate and it is the one of interest for the present analysis.
The figure shows that the first branch is a good approximation to waves in the continuous plate
if the element length is about one-tenth of the wavelength, that is approximately µx < 2π/10.
The second and third passing bands are due to the discrete FE model of the plate. It can be
noted that in the first passing band the model behaves as a low-pass filter for all the frequencies
below the cut-off frequency Ωa . The frequencies outside the passing bands, i.e. Ωa < Ω < Ωb
and Ω > Ωc , do not propagate. Although the presence of the bounding frequencies, Ωa , Ωb and
Ωc , is due to the discretisation of the continuous model, their evaluation can be of importance
since they are characteristics of the structure. In [13], for example, it has been proved that the
bounding frequencies for a 1-dimensional periodic system of symmetric elements, correspond
to the characteristic frequencies of the single element with its ends either fixed or free. Other
14
0.6
1500 Hz
800 Hz
0.4
200 Hz
µy
0.2
0
−0.2
−0.4
−0.6
−0.6
−0.4
−0.2
0
µx
0.2
0.4
0.6
Figure 5: Isotropic plate. Dispersion curves in the (µx , µy ) plane. Square element.
considerations about the bounding frequencies can be found in [11], [12], [14].
Figures 5 and 6 show the dispersion curves in the (µx , µy ) plane for the first propagating
band at different values of the propagation frequency. Two examples are considered. The first
concerns a square plate element of length L, Figure 5, while the second applies to a rectangular
element for which the distances between the nodes along x and y are taken as Lx = L and
Ly = 1/2L respectively, Figure 6. Since the plate is isotropic the wavenumber is independent
of the propagation direction and, as expected, the contour curves for the square plate element
are circular, Figure 5. On the other hand, the curves in Figure 6 relative to the rectangular
element are ellipses since µy = 1/2ky L and µx = kx L.
Some examples are now given in order to show the frequencies, f, and the associated wave
forms, V. These result from the eigenproblem defined in equation (11). The waveforms are
evaluated for plane waves that propagate in the directions θ = 0 and θ = π/4.
Illustrative example 1). Direction of propagation θ = 0, wavenumber 57.86 m−1 .
From equation 11, the frequencies and the corresponding eigenvectors are:
−0.0173
0
0
400
0
0
−0.001 + 0.1i .
(24)
f = 2.328 · 105 ; V =
5
−i
−0.001 + 0.1i
0
2.401 · 10
Illustrative example 2). Direction of propagation θ = π/4, wavenumber 57.86 m−1 .
15
0.4
1500 Hz
0.3
800 Hz
0.2
200 Hz
µy
0.1
0
−0.1
−0.2
−0.3
−0.4
−0.6
−0.4
−0.2
0
µx
0.2
0.4
0.6
Figure 6: Isotropic plate. Dispersion curves in the (µx , µy ) plane. Rectangular element.
From equation 11, the frequencies and the corresponding eigenvectors are:
−0.0244 0
0
400
i
−1 −1 .
f = 2.364 · 105 ; V =
5
−i
−1 1
2.364 · 10
(25)
The eigenvectors, V , correspond to the nodal displacements, which are translation in the z
direction, w, and rotations about the nodal x and y axes, θx = ∂w/∂y and θy = −∂w/∂x.
Since harmonic free waves propagate according to equation (21), it can be verified that, in the
two examples illustrated, only the first eigenvector represents a free propagating wave in the
continuous plate. The other two modes are due to the discretisation of the structure and do
not correspond to physical motion in the continuous plate.
Short movies files (available on request) have been produced with MATLAB in order to display
the wave form in the plate FE. The visualisation of the wave motion in the movies is carried
out for one square FE a cell of 6 × 6 square FEs. The movies show clearly that the first wave
forms correspond to a propagating disturbance for which the motion of every particle in the
plane perpendicular to the direction of propagation is the same.
3.1.2. Accuracy
Approximate numerical methods, such FEM, have been developed to numerically approximate
partial differential equations when is not possible to obtain an analytical solution. Nevertheless,
the simplifications and approximations used to obtain the model affect the output. One of the
FE approximation errors consists in dispersion errors that are caused by the difference between
the numerical evaluation of the wavenumber and the true wavenumber. In order to evaluate the
16
+
Re(µ ) Numerical
+
Im(µ ) Numerical
Analytical
7/4
3/2
5/4
µ+/π
1
3/4
1/2
1/4
0
0
2
4
6
8
10
Ω
12
14
16
18
20
Figure 7: Isotropic plate. Complex–valued dispersion curve for θ = 0.
dispersion error, as an example, a comparison of the analytical and the numerical dispersion
curves is given in Figure 7. The figure shows the nondimensional wavenumber µ+ = kL against
Ω for a wave propagating in the positive x direction. Figure 7 shows that the numerical solution
is close to the analytical one until µ = π. In term of wavelength this means L < λ/2. The
same considerations apply to any value of the propagation direction θ, that is Lx < λx /2 and
Ly < λy /2. Hence, it can be concluded that the accuracy of the method is good when common
requirements for FE discretisation are satisfied: Lx < λx /10 and Ly < λy /10. This is an oftenused criterion that can be used to ensure that there are enough elements per wavelength to
accurately characterise the wave.
The dispersion errors also depend on the propagation direction. Figure 8 shows the relative
error in the estimated frequency and in the estimated wavenumber as a function of θ, Figure
8.(a) and Figure 8.(b) respectively.
Different cases are considered: a square
√ plate of length Lx = L and rectangular plates for
which Lx = L, Ly = 1/2L and Ly = 2L. As an example, the frequencies are obtained for
a plane wave whose nondimensional wavenumber is kL = 0.2893 while the wavenumbers are
obtained for nondimensional frequency Ω = 0.0837.pAs expected, the errors in the frequency
and wavenumber evaluation increase as the distance L2x + L2y between the nodes 1−4 increases
and reach their maximum as the propagation direction approaches π/4 + mπ/2, m integers.
The fact that the plot in Figure 8.(b) is piecewise constant is due to the methods, (Newton
method and Powell’s method ), used to solve the transcendental equation in section 2.5. In fact,
when these methods are applied to find a root of complex functions, the set of points where
the methods converge to that root is a domain of attraction.
17
0.7
L =1/2L
y
L =L
y 1/2
L =2 L
(a)
y
0.6
error %
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
θ/π
0.25
Ly=1/2L
L =L
y 1/2
Ly=2 L
(b)
0.2
error %
0.15
0.1
0.05
0
0
0.2
0.4
0.6
0.8
1
θ/π
1.2
1.4
1.6
1.8
2
Figure 8: Isotropic plate. Relative error in the estimated frequency, figure (a), and in the wavenumber,
figure (b), as a function of the wave propagation direction θ, kL = 0.2893, Ω = 0.0837.
3.1.3. Complex–valued propagation constants
In this section some examples are given in order to show the solutions obtained from equations
(15) and (17) for a given frequency and a given propagation direction. A square plate element
of length L is considered. As an example, the solutions are evaluated at Ω = 0.0837 for different
values of θ. Figure 9 shows the location of the nondimensional wavenumber µ = kL obtained
−iσ
solving equation (15) for θ = arctan(1/2).
q Since in equation (15), γ = e , the nondimensional
wavenumbers are evaluated as kL = (µ2x + µ2y ), where µx = m1 σ, µy = m2 σ and σ = i ln(γ).
As another example, Figure 10 shows the nondimensional wavenumbers obtained from problem
(17) for θ = π/3.
Equation (15) can also be arranged to obtain µy for a given µx (or µx for a given µy ) for a given
18
Numerical
Analytical
4
3
2
Im(kL)
1
0
−1
−2
−3
−4
−8
−6
−4
−2
0
2
4
6
8
Re(kL)
Figure 9: Isotropic plate. Location of the nondimensional wavenumber –obtained solving equation (15)–
in the complex (Re(kL), Im(kL)) plane. θ = arctan(1/2), Ω = 0.0837.
frequency as a function of θ. Figure 11 shows the real and imaginary part of the solutions µy ,
for given µx as a function of the propagation direction. Solving equation (15), 6 solutions are
obtained which corresponds to 3 pairs of waves which either decay or propagate in the negative
and positive y direction. Since damping is not considered, the solutions for freely propagating
waves satisfy the equation |e−iµy | = 1, hence solutions 1 and 2 represent the component along
y of the wave that propagates along ±θ direction. Solutions 3, 4, 5 and 6 represent the y
components of evanescent waves with amplitudes that decrease in the θ and −θ directions
respectively. The real part of µy shows that adjacent nodes along the y direction vibrate in
phase for solution 3 and 4 and in counter phase for solution 5 and 6.
3.1.4. Sensitivity analysis of eigenvalues
From the analysis carried out, it is clear that there is a set of solutions which do not correspond
to the results for the continuous structure but are due to the discretisation. Thus there is need
for a systematic procedure that enables results which represent the behaviour of the continuous
structure to be evaluated. This means that the eigenmodes that correspond to waves in the
continuous model must be identified. A possible criterion is to look at the sensitivity of the
results to a particular parameter. Since a design parameter in the FE model is the distance
between the nodes, a small change in this parameter should result in relatively small changes
for the results representative of wave motion in the continuous structure and in relatively large
changes in those associated with the discrete structure.
A direct approach to obtain the eigenvalue sensitivities is to differentiate equation (11) with
respect to the element length L. If the eigenvalues are simple (not multiple) and the eigenvectors
19
Numerical
Analytical
4
3
2
Im(kL)
1
0
−1
−2
−3
−4
−10
−8
−6
−4
−2
0
2
4
6
8
10
Re(kL)
Figure 10: Isotropic plate. Location of the nondimensional wavenumber –obtained solving equation
(17)– in the complex (Re(kL), Im(kL)) plane. θ = π/3, Ω = 0.0837.
are normalised with respect to the mass matrix, the expression for the derivative of equation
(11) is
∂K
∂ωi2
T
2 ∂M
qi ;
i = 1, . . . , n.
(26)
= qi
− ωi
∂L
∂L
∂L
In simple cases the stiffness and mass matrices can be evaluated as functions of the dimension
of the element and their derivatives with respect to L can be determined analytically. If only
a numerical evaluation of the matrices is available, the derivative of the element matrices may
be approximated by its first order numerical derivative as
∂K(L)
∂M (L)
K(L + ∆L) − K(L)
M (L + ∆L) − M (L)
=
;
=
.
(27)
∂L
∆L
∂L
∆L
For example, for a square plate element of dimension L = 0.005 m, the nondimensional analytical variation of the frequency with respect to the L is:
∂Ω22
∂Ω23
∂Ω21
= 5.433 · 10−5 ;
= −1.954 · 106 ;
= −2.018 · 106 .
(28)
∂L
∂L
∂L
From (28) shows clearly that the second and third modes are very sensitive to a change in the
element dimension and therefore are due to the discretisation process.
3.1.5. Sensitivity analysis of the wavenumbers
The FE estimates depend on the element aspect ratio. This is particularly significative for the
higher periodic structure branches while this is not true for solutions which represent wave–
modes in the continuous model. Figures 12 and 13 show, as an example, the wavenumbers for
20
5/4
(a)
5
6
1
y
Re(µ )/π
3/4
1/2
1/4
2
1
3
4
0
−1
−3/4
−1/2
−1/4
0
1/4
θ/π
1/2
3/4
1
3/4
1
(b)
1/2
6
1/4
y
Im(µ )/π
4
1 2
0
3
−1/4
5
−1/2
−1
−3/4
−1/2
−1/4
0
θ/π
1/4
1/2
Figure 11: Isotropic plate. Real part, figure (a), and imaginary part, figure (b), of µy for given µx and
given frequency as a function of θ.
various element aspect ratios. The complex values of the nondimensional wavenumbers, kL,
are obtained solving the eigenproblem (17) for θ = π/3 and Ω = 0.0837.
Figure 12 shows that while solutions 1, 2, 3 and 4 seem to be not affected by small changes in
the element geometrical parameters, the other solutions are very sensitive to changes in the
aspect ratio.
Figure 13 shows the sensitivity in the estimated wavenumbers as a function of L. The error
in the wavenumber is evaluated as (k(L0 ) − k)/k(L0 )% where L0 is a given initial length. The
figure shows that solutions 1, 2, 3 and 4 are insensitive to the variation of L while a change of
1% in the length of the element results in a change of nearly 1% for the other solutions. It can
be concluded that solutions 1, 2, 3 and 4 can be considered as the solutions representative of
the wave motion in the continuous model.
21
5
L =1/2L
y
L =L
y 1/2
L =2 L
4
y
3
2
Im(kLx)
1
4
1
3
0
2
−1
−2
−3
−4
−5
−8
−6
−4
−2
0
2
4
6
8
Re(kLx)
Figure 12: Isotropic plate. Location of the solutions in the complex (Re(kLx ), Im(kLx )) plane,
Ω = 0.0837, θ = π/3.
Solution 1, 2, 3, 4
Other solutions
0.5
0
[k(L )−k(L)]/k(L )%
1
0
0
−0.5
−1
0.99
1
1.01
L/L0
Figure 13: Isotropic plate. Sensitivity analysis of the wavenumber with respect to change in plate
dimension.
22
3.2. Orthotropic plate
In this section an orthotropic plate made of a Glass/Epoxy composite material with the fibres
oriented in the x direction is considered. The plate has uniform thickness and the material
properties are as follows: Young’s modulus in the fibre direction Ex = 39GP a, Young’s modulus in the transverse direction Ey = 8.6GP a, in-plane shear modulus Gxy = 3.8GP a, major
Poisson’s ratio ν = 0.28, density ρ = 2100kg/m3 .
For such a plate the equation for bending vibrations is
Dx
∂4w
∂4w
∂2w
∂4w
+
2H
+
D
=
−ρh
,
y
∂x4
∂x2 ∂y 2
∂y 4
∂t2
(29)
where h is the plate thickness and H = D1 + 2Dxy . The rigidities are defined as
Dx =
Ex h3
;
12(1 − ννm )
Dy =
Ey h3
;
12(1 − ννm )
D1 =
Ey νh3
;
12(1 − ννm )
Dxy =
Gxy h3
,
12
where νm is the minor Poisson’s ratio. For more details see [29]. Substituting (21) into (29)
and considering that kx = k cos θ and ky = k sin θ, the free wavenumber is
s
√ 4
ρh
k= ω
.
(30)
4
Dx cos θ + 2H cos θ2 sin θ2 + Dy sin θ4
A nondimensional frequency is defined as
Ω=
ωL2x
q
e
ρh/D,
e = Ex h3 /12(1 − ννm ). As in the previous section, if nothing different is specified the
where D
plate element is assumed to be square and of length L.
3.2.1. Real–valued dispersion curves
Figure 14 shows the surface plot of the nondimensional frequencies with respect to µx = kx L
and µy = ky L. As expected the surfaces are symmetric with respect to the µx = 0 and µy = 0
axes.
Figures 15 shows the dispersion curves for the first propagating band in the (µx , µy ) plane
at different values of the frequency. Although the element is square, the contour curves are
ellipses because the material is orthotropic. Figure 16 shows the nondimensional wavenumber
µ+ = kL against Ω for a wave propagating in the positive x direction. The Figure shows that
the results are accurate if the length of the element is small compared with the wavelength.
The same conclusion applies to any value of the propagation direction θ, that is Lx < λx /10
and Ly < λy /10.
23
10
8
Ω1
6
4
2
0
1
1
0.5
0.5
0
0
−0.5
−0.5
−1
µ /π
−1
µx/π
y
25
Ω2
20
15
10
5
1
1
0.5
0.5
0
0
−0.5
−0.5
−1
µ /π
−1
µ /π
y
x
55
50
45
40
Ω
3
35
30
25
20
15
10
1
1
0.5
0.5
0
0
−0.5
µ /π
−0.5
−1
−1
y
µ /π
x
Figure 14: Orthotropic plate. Propagation surfaces.
24
1
1500 Hz
0.8
800 Hz
0.6
200 Hz
0.4
µy
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
0.6
0.4
0.2
0
−0.2
µ
−0.4
−0.6
x
Figure 15: Orthotropic plate. Dispersion curves in the (µx , µy ). Square element.
+
Im(µ ) Numerical
+
Re(µ ) Numerical
Analytical
7/4
3/2
5/4
µ+/π
1
3/4
1/2
1/4
0
0
2
4
6
8
10
Ω
12
14
16
18
Figure 16: Orthotropic plate. Complex–valued dispersion curve for θ = 0.
25
20
3.2.2. Complex–valued propagation constants
In this section an example is given in order to show the solutions obtained from equation (17)
for a given frequency and a given propagation direction. Figure 17 shows the location of the
nondimensional wavenumber in the complex (Re(kL), Im(kL) plane for several values of θ.
The solutions are evaluated for Ω = 0.1003. It can be concluded that solutions 1, 2, 3 and 4 are
the solutions representative of the wave motion in the continuous model while the others are
artifacts of the discretisation of the structure.
θ=π/3 Numerical
θ=π/4 Numerical
θ=π/6 Numerical
θ=π/6, π/4, π/3 Analytical
4
3
2
1
Im(kL)
4
3
0
1
2
−1
−2
−3
−4
−8
−6
−4
−2
0
2
4
6
8
Re(kL)
Figure 17: Orthotropic plate.
Location of the nondimensional wavenumber in the complex
(Re(kL), Im(kL)) plane, Ω = 0.1003.
26
4. Wave propagation in curved and cylindrical shells
Thin cylindrical shells are used in many engineering applications and their modal characteristics
have been a subject of study of great interest. The free vibration characteristics of isotropic
cylindrical shells have been obtained by various approximate theories. A good summary of
these theories is given by Leissa in [30]. There have been also a number of previous studies
about wave propagation in thin cylindrical shells and a brief summary of some of them is given
here.
D. C. Gazis analysed the propagation of free harmonic waves in a thin long cylinder. He
provided the frequency equation for an arbitrary number of waves around the circumference
in [31] and he presented the results for free harmonic waves in [32]. Kumar and Sthepens,
[33], studied dispersion of flexural waves in isotropic shells of various wall thickness using the
exact three dimensional equations of linear elasticity. In [34], Fuller considered the effect of wall
discontinuities on flexural waves in an isotropic cylindrical shell. He also provided the dispersion
curves in a semi-infinite thin walled shell for circumferential modal number n = 1. In [35], Fuller
and Fahy analysed the dispersion and energy distributions of free waves in isotropic cylindrical
shells filled with fluid. The dispersion curves have been obtained for water-filled steel and hard
rubber shells vibrating in the circumferential modes of order n = 0 and n = 1. More recently
Langley, [36, 37], extended the analysis to general helical motion in isotropic cylindrical shells
and studied the modal density and mode count of the out-of-plane modes for isotropic thin
cylinders and curved panels. Tyutekin, [38], obtained the dispersion curves in an isotropic
cylindrical shell for different angles of propagation of the helical waves. He also evaluated the
dispersion characteristics and the eigenfunctions of circumferential waves for the Neumann and
Dirichlet boundary conditions, [39].
w3
α
w4
z
η4
4
η3
ξ4
y
w2
α
w1
η1
3
x
1
ξ3
η2
ξ1
ξ2
h
2
R
α
α
α
(a)
(b)
Figure 18: Segment of cylinder, figure (a), and rectangular flat shell elements, figure (b).
The aim of this section is to apply the technique proposed in the previous sections to the case of
cylindrical shells. An isotropic cylindrical shell is considered in the first part. The real–valued
dispersion curves are obtained for different circumferential mode numbers and considerations
about the nature of the propagating waves are provided. In the second part, the same results
are extended to the more complicated case of a sandwich cylindrical shell.
27
In the following analysis the FE-software ANSYS is used to obtain the mass and stiffness
matrices of the shell finite element. The element SHELL63 is used for the case of the isotropic
shell while the element SHELL181 is used for the more general case of a sandwich shell. The
elements are rectangular four node shell elements. They have both bending and membrane
capabilities and they have six degrees of freedom at each node: translations in the nodal x, y,
and z directions and rotations about the nodal x, y, and z axes. The nodes are numbered as in
Figure 18. These elements are flat, hence the curved surface is approximated by piece–wise–flat
surfaces and thus, adjacent elements, as shown in Figure 18.(b), are identical flat elements
except the normal is rotated through the angle α. The accuracy in modelling the shells is
governed by the Mindlin-Reisner shell theory. Figure 18 shows a segment of the cylinder and
the rectangular flat shell elements: R represents the radius of curvature and h is the shell
thickness; x, y, z represent the global coordinates in the circumferential, the axial and the
normal directions respectively while ξ, η, w represent the nodal coordinates.
Harmonic free waves in the FE cylindrical shell that propagate over an element have the form
φ = Aei(ωt−µx −µy ) ,
(31)
where µy , µx are the propagation constants along the axial and circumferential directions.
Equation (31) allows any value of µx and it represents a disturbance that propagate in a
helical pattern. It can be noted that equation (31) is the same used for the plate case. This
may be interpreted by considering that the wave–FE model works for the curved shell as an
equivalent “flat model” when the flat shell is infinite in both x and y directions. If k is the
helical wavenumber and θ is the wave propagation direction with respect to the x axis, the
wavenumber can be represented by its components kx = k cos(θ) and ky = k sin(θ) in the x and
y directions. Hence µx = kx Lx and µy = ky Ly where Lx and Ly are the element dimensions.
For an uniform cylindrical shell, due to the closure of the cylindrical geometry, periodic values
in the angle α = x/R are used and the solution of the partial differential equations of motion
is evaluated for integer multiples n of the angle α. The values of the number n indicates what
in the literature is referred to as circumferential modes. Since µx identifies the change of phase
between one node and the next one along the circumferential direction, to obtain a specific
value of µx at all nodes for the cylindrical geometry, then N µx = n2π, where N = 2πR/Lx is
the number of FEs in the circumferential direction. Therefore µx = nLx /R is the propagation
constant in the x direction that corresponds to the circumferential mode number n.
The shell elements in ANSYS are flat elements. In order to model the desired curvature, the
DOFs of the nodes 2 and 3 are rotated, with respect to the global coordinates x, y, z, through an
angle α to get DOFs of the nodes 1 and 4 of the next element, Figure 18.(b). Then periodicity
conditions are applied. Therefore the nodal degrees of freedom are transformed by the matrix
I 0 0 0
0 r 0 0
R=
(32)
0 0 r 0 ,
0 0 0 I
where 0 and I are zero and identity matrices of order 6 × 6 while r refers to the affine transfor-
28
mation
r=
cos(α)
0
sin(α)
0
0
0
0 − sin(α)
0
0
0
1
0
0
0
0
0 cos(α)
0
0
0
0
0
cos(α) 0 − sin(α)
0
0
0
1
0
0
0
sin(α) 0 cos(α)
.
(33)
The new mass and stiffness matrices with respect to the global system of coordinates become
f = RT MR and K
e = RT KR. This introduces coupling between the out–of–plane and the
M
in–plane behaviour of the shell.
An important value in the wave behaviour of a thin cylindrical shell is the ring frequency.
The ring frequency is the frequency at which a wave travelling at the longitudinal propagation
speed in a plate has a wavelength equal to the circumference. This frequency divides the shell
behaviour into two regions: above the ring frequency the shell behaves as a flat plate while near
and below the ring frequency the the effect of the curvature stiffen the structure and results in
a more complicated behaviour.
4.1. Isotropic shell
In the following results the nondimensional frequency is defined as Ω = ω/ωr where
p
ωr = 1/R E/ρ(1 − ν 2 )
is the ring frequency for the cylindrical shell. The frequency range of interest in the present
analysis includes the ring frequency.
A steel cylindrical shell is considered. The nondimensional thickness of the shell is h = h/R =
0.05 while the material properties are: Young’s modulus E = 19.2 · 1010 N/m2 , Poisson’s ratio
ν = 0.3, density ρ = 7800 kg/m3 .
Figures 19–20 show the real-valued dispersion curves for the circumferential modes of order
n = 0, 1, 2, and 3. Branches 1, 2, 3 are indicated. Since there are 6 DOFs per node, six real
passing bands are obtained. However, only the first three are shown since they are the ones that
correspond to propagating waves in the continuous shell model while the other are numerical
artefacts due to the discretisation.
These wave modes can be classified as primarily longitudinal, torsional or flexural. The longitudinal and torsional waves correspond to non–dispersive waves with phase speeds in the flat
plate
s
s
E
G
and cT =
,
cL =
2
ρ(1 − ν )
ρ
respectively, while the flexural waves are dispersive with a wave speed
s
√ 4
Eh3
cB = w
.
12ρ(1 − ν 2 )
Branches 1, 2 and 3 broadly correspond to flat–plate flexural, torsional and extensional waves
respectively. It can be noted that the minimum cut–off frequency for branches 1, 2, and 3 equal
29
0.08
0.07
1
0.06
1
µ
y
0.05
0.04
0.03
2
0.02
3
2
3
0.01
0
0
0.2
0.4
0.6
0.8
1
1.2
Ω
1.4
1.6
1.8
2
Figure 19: Isotropic shell: real–valued dispersion curves. · · · n = 0, — n = 1.
0.08
0.07
1
0.06
1
µy
0.05
0.04
0.03
2
2
0.02
3
3
0.01
0
0
0.5
1
1.5
Ω
2
2.5
3
3.5
Figure 20: Isotropic shell: real–valued dispersion curves. · · · n = 2, — n = 3.
30
0.2
branch 1
branch 2
branch 3
0.18
0.16
0.14
µy
0.12
flexural
0.1
0.08
torsional
0.06
longitudinal
0.04
0.02
0
0
0.5
1
1.5
2
2.5
Ω
3
3.5
4
4.5
5
Figure 21: Isotropic shell: dispersion curves for n = 0. Continuous lines refer to extensional, torsional
and flexural wave speed in the flat shell.
the frequency at which a wave travelling at the flexural, torsional and extensional propagation
speed in a plate has a wavelength equal to the circumference. In fact, the n = 0 branch 3
cuts-off at the frequency Ω = 1, the n = 1 branch 2 cuts-off at Ω = 0.6 and the n = 2 branch
1 cuts–off at the frequency Ω = 0.04, see Figures 19 and 20. As an example, Figure 21 shows
a graphical comparison between the numerical results for branches 1, 2 and 3 (n = 0) and the
analytical dispersion curves for the flexural, torsional and extensional waves in a thin flat plate.
It can be seen that above the ring frequency the shell behaves as a flat plate. The difference
between these branches and the flexural, torsional and extensional dispersion curves becomes
smaller as the circumferential number increases.
Figures 22–25 show the wave speed of the branches 1, 2 and 3 where the nondimensional wave
speed c∗ is defined as the ratio between the wave speed and the extensional wave speed cL . For
n = 0 only waves 1 and 2 exist for small value of the axial propagation constant, that is below
the ring frequency. They propagate as in–plane, extensional and torsional waves respectively.
Hence, for Ω < 1 the cylindrical shell behaves dynamically as a membrane. For n = 1, both
waves 1 and 2 propagate for small value of µy . However, their behaviour cannot be considered
as purely flexural, torsional or extensional. For increasing values of the circumferential wave
number, below the ring frequency, only wave 1 exists and its behaviour approaches that of a
flexural wave in an infinite plate.
The dispersion curves can also be represented in the (µx , µy ) plane for fixed values of Ω. As
examples, Figures 26–28 show the curves representing dispersion branches for Ω = 0.1, 0.2,
Ω = 0.5, 1 and Ω = 1.5. It can be seen in Figures 26–27 that only waves 1 and 2 propagate at
these frequencies while wave 3 cuts-on at Ω = 1. In particular for wave 1, for low values of the
frequency, there exist regions in which a value of µy corresponds to two distinct values of µx .
31
2
branch 1
branch 2
branch 3
1.8
n=0
1.6
1.4
1.2
c
*
extensional wave speed
1
0.8
torsional wave speed
0.6
0.4
flexural wave speed
0.2
0
0
0.02
0.04
0.06
µy
0.08
0.1
0.12
Figure 22: Isotropic shell: nondimensional wave speed for n = 0. Continuous lines refer to extensional,
torsional and flexural wave speed in the flat shell.
2
branch 1
branch 2
branch 3
1.8
n=1
1.6
1.4
1.2
c
*
extensional wave speed
1
0.8
torsional wave speed
0.6
0.4
flexural wave speed
0.2
0
0
0.02
0.04
0.06
µ
0.08
0.1
0.12
y
Figure 23: Isotropic shell: nondimensional wave speed for n = 1. Solid lines refer to extensional,
torsional and flexural wave speed in the flat shell.
32
1.8
branch 1
branch 2
branch 3
1.6
n=2
1.4
1.2
extensional wave speed
c
*
1
0.8
torsional wave speed
0.6
0.4
flexural wave speed
0.2
0
0
0.02
0.04
0.06
0.08
0.1
0.12
µy
Figure 24: Isotropic shell: nondimensional wave speed for n = 2. Solid lines refer to extensional,
torsional and flexural wave speed in the flat shell.
1.8
branch 1
branch 2
branch 3
n=3
1.6
1.4
1.2
extensional wave speed
c
*
1
0.8
torsional wave speed
0.6
0.4
flexural wave speed
0.2
0
0
0.02
0.04
0.06
µ
0.08
0.1
0.12
y
Figure 25: Isotropic shell: nondimensional wave speed for n = 3. Solid lines refer to extensional,
torsional and flexural wave speed in the flat shell.
33
0.02
0.015
1
0.01
1
y
0.005
2
µ
0
−0.005
−0.01
−0.015
−0.02
−0.05
−0.04
−0.03
−0.02
−0.01
0
µx
0.01
0.02
0.03
0.04
0.05
Figure 26: Isotropic shell: dispersion curves in the (µx , µy ) plane. — Ω = 0.1, · · · Ω = 0.2.
0.06
0.04
1
1
µy
0.02
2
0
2
3
−0.02
A
B
−0.04
−0.06
−0.08
−0.06
−0.04
−0.02
0
µx
0.02
0.04
0.06
0.08
Figure 27: Isotropic shell: dispersion curves in the (µx , µy ) plane. — Ω = 0.5, · · · Ω = 1.
34
0.1
0.08
1
0.06
0.04
2
µy
0.02
3
0
−0.02
−0.04
−0.06
−0.08
−0.1
−0.1
−0.05
0
0.05
µx
0.1
Figure 28: Isotropic shell: dispersion curves in the (µx , µy ) plane for Ω = 1.5.
As an example, points A and B are shown in Figure 27 for a given value of µy . A discussion
about the nature of the two different waves associated with the lower and higher values of µx
for given value of µy can be found in [36]. A graphical interpretation of the energy flow for
points A and B can be made using the dispersion curves. The group velocity is defined by
cg =
dω
.
dk
(34)
Since the dispersion curves in the (µx , µy ) plane have been obtained as contour curves such
that Ω = f (µx , µy ), the gradient is always normal to the curves Ω = constant. The gradient
∇Ω is a vector whose components are the partial derivatives of Ω, that is
∂Ω ∂Ω
,
∇Ω =
µx µ y
T
.
(35)
Therefore, considering (34) and (35), it can be concluded that the direction of the group velocity
(or the direction of the energy flow) is the normal to the contours. It can be seen in Figure 27
that the normal vector to the curve at point B has a negative component with respect to the
x direction and therefore it has a negative group velocity in the circumferential direction.
For Ω = 1.5, Figure 28, all the three waves propagates with a dynamic behaviour that is similar
to that of the extensional, torsional and flexural waves in an infinite plate. It is of note that all
these curves reveal the “anisotropic” properties of the curved shell element with respect to the
corresponding flat shell element.
35
4.2. Sandwich shell
Figure 29: Sandwich shell: layers stacking and material properties.
The sandwich shell is made of a bottom skin laminate, a foam core and a top skin laminate.
The upper and lower laminae are constituted by 4 orthotropic sheets of glass/epoxy. Each
skin has a lay–up of [+45/-45/-45/+45] and a total thickness of 4 mm. The core is a 10 mm
polymethacrylamide ROHACELL foam with 110WF of density. The different layers and the
material properties are shown in Figure 29. The nondimensional thickness is h = 0.018.
Real–valued dispersion curves are shown in Figures 30–31 for the circumferential modes n =
0, 1, 2 and 3. The ring frequency can be considered as the first transition frequency for the thin
cylinder and it may be evaluated as the frequency at which the third wave starts propagating.
Herein, since for n = 0 the branch 3 cuts-on at 617Hz, this frequency value is taken as the ring
frequency for the sandwich cylindrical shell. The minimum cut-off frequencies for branches 2
and 1 are 550.3Hz and 11.5Hz respectively. As expected the wave behaviour below the ring
frequency is very complex and cannot be described simply in terms of torsional, extensional
and flexural waves alone. It can be seen in Figure 30 that the n = 1 branches 1 and 2 exhibit
more than one possible value of µy for the same value of frequency. For the n = 1 branch 1,
three values of µy can be identified by the same frequency when 384Hz< f < 399Hz, Figure
30. The lower and the higher value correspond to waves with positive group velocity in the y
direction while the middle value corresponds to a wave which has a negative group velocity in
the y direction. For n = 1 branch 2, in the frequency range 428Hz< f < 550Hz, two values of
µy correspond to two different waves travelling in opposite y direction, Figure 30. In particular,
the wave associated with the lower value of µy has a negative group velocity in the y direction.
Similar conclusions can be drawn for the dispersion curves in Figure 31. For n = 0 and n = 1
only waves 1 and 2 propagate below the ring frequency. For n = 2, wave 2 still propagates for
frequencies close to the ring frequency, while for n = 3 only wave 1 propagates below the ring
frequency.
36
0.06
0.05
µy
0.04
0.03
1
0.02
2
0.01
2
3
0
0
200
400
600
800
3
1000
1200
1400
f [Hz]
Figure 30: Sandwich shell: real–valued dispersion curves. · · · n = 0, — n = 1.
0.1
0.09
0.08
0.07
µy
0.06
0.05
1
0.04
0.03
2
0.02
2
0.01
3
0
0
200
400
600
800
1000
1200
1400
1600
f [Hz]
Figure 31: Sandwich shell: real–valued dispersion curves. · · · n = 2, — n = 3.
37
0.06
0.04
1
0.02
1
µ
y
2
2
0
−0.02
−0.04
−0.06
−0.08
−0.06
−0.04
−0.02
0
µx
0.02
0.04
0.06
0.08
Figure 32: Sandwich shell: dispersion curves in the (µx , µy ) plane. — f = 200Hz, · · · f = 500Hz.
0.1
0.08
1
0.06
0.04
µ
y
0.02
2
0
3
−0.02
−0.04
−0.06
−0.08
−0.1
−0.1
−0.05
0
µ
0.05
0.1
x
Figure 33: Sandwich shell: dispersion curves in the (µx , µy ) plane for f = 800Hz.
38
Low frequency dispersion curves in the (µx , µy ) plane are shown in Figures 32–34 for different
values of frequency. It can be seen in Figures 32–34 that there exists regions for branch 2
in which distinct values of µx correspond to the same value of µy and distinct values of µy
correspond to the same value of µx . For every one of these points, considerations about the
energy flow can be made analogous to that for the isotropic shell. A comparison between the
dispersion curves 32–34 and results obtained by Heron using an analytical approach for a similar
curved panel [40] has shown that the dispersion curves are almost identical in shape.
0.15
1
0.1
0.05
2
µ
y
3
0
−0.05
−0.1
−0.15
−0.1
−0.05
0
µx
0.05
0.1
0.15
Figure 34: Sandwich shell: dispersion curves in the (µx , µy ) plane for f = 1500Hz.
39
5. Concluding remarks
A method for the numerical prediction of the wave characteristics of 2–dimensional structures
using standard FEM has been discussed. With this method only one rectangular FE with 4
nodes is used to describe the dynamics of the system reducing drastically the cost of calculations.
The mass and stiffness matrices of the element are found using conventional FE approaches and
thus the output of commercial FE package can be used. The mass and stiffness matrices are
subsequently post–processed using periodicity conditions in order to formulate an eigenproblem whose eigenvalues and eigenvectors describe the wave propagation characteristics. This
eigenproblem is the kernel of the method and can be cast in a form to provide the propagating
frequencies and the waveforms for given values of the propagation constants, or to evaluate the
propagation constants for given values of the frequency. The latter approach predicts complex
wavenumbers and therefore can be extended to the analysis of general cases, for example in
presence of energy dissipation.
The method has been applied to several examples. The first example concerns the out–of–
plane vibration of an isotropic plate. Using this example some considerations about general
features of the method have been drawn. It has been found that the frequency of propagation
is a periodic function of the real propagation constants. It has also been observed that the
waves propagate across the structure only within some propagating bands. The number of
these bands, and therefore the number of frequencies and waveforms, depend on the number
of the nodal DOFs chosen for the element. Further, it has been shown that some of them
represent the behaviour of the continuous model, while others are numerical solutions due to
the discretisation of the system. A sensitivity analysis has been carried out in order to evaluate
the set of solutions which correspond to the continuous structure. The second example deals
with an orthotropic plate made of Glass/Epoxy composite material. In the last part cylindrical
shells have been considered. In particular an isotropic and a sandwich curved shell have been
analysed. These last two cases have been studied by postprocessing an ANSYS FE model.
Since the shell elements obtained from ANSYS are flat shell elements, the matrices of a single
element are manipulated to model the desired curvature for the corresponding curved shell. The
example of the sandwich cylindrical shell has proved the usefulness of the method for predicting
dispersion curves when analytical solutions are not available.
The main advantages of the technique can be thus summarised as follows:
• the computational cost becomes independent of the size of the structure since the method
involves typically just one finite element to which periodicity conditions are applied;
• the technique is very flexible since standard FE routines can be used and therefore a wide
range of structural configurations can be easily analysed;
• finally, the propagation constants for plane harmonic waves can be predicted for different
propagation directions along the structure.
The method is seen to give accurate predictions when the length of the element is small enough
compared to the wavelength, that is when the length of the single element is less then 1/6–1/10
of the wavelength.
Further work should be done to improve the analysis presented in this memorandum. In particular, advances of the post–processing algorithms are needed in order to obtain the complete set
of the complex–valued dispersion curves. Future work may also concern the evaluation of the
40
response to external loads such as point forces or random pressure fields. It is also of interest
to extend the applicability of the method to non–rectangular finite elements.
41
References
[1] S. Finnveden, “Exact spectral finite element analysis of stationary vibrations in a railway
car structures”, Acta Acustica, Vol. 2, pp. 461-482, 1994.
[2] L. Gavric, “Computation of propagative waves in free rail using a finite element technique”,
Journal of Sound and Vibration, Vol. 184, pp. 531-543, 1995.
[3] U. Orrenius, S. Finnveden, “Calculation of wave propagation in rib-stiffened plate structures”, Journal of Sound and Vibration, Vol. 198, pp. 203-224, 1996.
[4] S. Finnveden, “Spectral finite element analysis of the vibration of straight fluid-filled pipes
with flanges”, Journal of Sound and Vibration, Vol. 199, pp. 125-154, 1997.
[5] P. J. Shorter, “Wave propagation and damping in linear viscoelastic laminates”, Journal of
the Acoustical Society of America, Vol. 115, pp. 1917-1925, 2004.
[6] J. F. Doyle, “Wave propagation in structures”, Springer-Verlag, 1997.
[7] J. R. Banerjee, “Dynamic stiffness formulation for structural element: a general approach”,
Computers Structures, Vol. 63(1), pp. 101-103, 1997.
[8] H. von. Flotow, “Disturbance propagation in structural networks”, Journal of Sound and
Vibration, Vol. 106, pp. 433-450, 1986.
[9] L. S. Beale, M. L. Accorsi, “Power flow in two- and three-dimensional frame strucutres”,
Journal of Sound and Vibration, Vol. 185, pp. 685-702, 1995.
[10] L. Brillouin, “Wave propagation in periodic structures”, Dover Publications, 1953.
[11] D. J. Mead, “Free wave propagation in periodically supported infinite beams”, Journal of
Sound and Vibration, Vol. 11(2), pp. 181-197, 1970.
[12] R. M. Orris, M. Petyt “A finite element study of harmonic wave propagation in periodic
structures”, Journal of Sound and Vibration, Vol. 33(2), pp. 223-236, 1974.
[13] D. J. Mead, “wave propagation and natural modes in periodic system: I. mono-coupled
systems”, Journal of Sound and Vibration, Vol. 40(1), pp. 1-18, 1975.
[14] D. J. Mead, “Wave propagation and natural modes in periodic systems: II. Multi-coupled
systems, with and without damping”, Journal of Sound and Vibration, Vol. 40(1), pp.
19-39, 1975.
[15] D. J. Mead, “Wave propagation in continuous periodic structures: research contributions
from Southampton”, Journal of Sound and Vibration, Vol. 190, pp. 495-524, 1996.
[16] D. Duhamel, B. R. Mace, M. J. Brennan, “Finite element analysis of the vibration of
waveguides and periodic structures”, Journal of Sound and Vibration, Vol. 294, pp. 205220, 2006.
[17] B. R. Mace, D. Duhamel, M. J. Brennan, L. Hinke, “Finite element prediction of wave
motion in structural waveguides”, Journal of the Acoustical Society of America, Vol. 117,
pp. 2835-2843, 2005.
42
[18] B. R. Mace, D. Duhamel, M. J. Brennan, L. Hinke, “Wavenumber prediction using finite element analysis”, Eleventh International Congress on Sound and Vibration, St. Petersburg,
2004.
[19] A. Y. A. Abdel-Rahman, “Matrix analysis of wave propagation in periodic systems”, PhD
thesis, University of Southampton, 1979.
[20] R. A. Tapia, “The Kantorovich Theorem for Newton’s Method”, The American Mathematical Monthly, Vol. 78(4), pp. 389-392, 1971.
[21] K. V. Singh, Y. M. Ram, “Transcendental eigenvalue problem and its applications”, AIAA
Journal, Vol.40(7), pp. 1402-1407, 2002.
[22] W. H. Yang, “A method for eigenvalues of sparse λ-matrices”, International Journal for
Numerical Methods in Engineering, Vol.19(6), pp. 943-948, 1983.
[23] M. J. D. Powell, “A Fortran Subroutine for Solving Systems of Nonlinear Algebraic Equations”, Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, Ch.7, 1970.
[24] L. E. Bateson, M. A. Kelamanson, C. Knudsen “Solution of a transcendental eigenvalue
problem via interval analysis”, Computers and Mathematics with Applications, Vol.38, pp.
133-142, 1999.
[25] R. E. Moore, “Interval Analysis”, Prentice-Hall, Englewood Cliff, N. J., 1966.
[26] E. Kreyszig, “Advanced Engineering Mathematics”, John Whiley Son, 1999.
[27] C. F. Gerald, “Applied Numerical Analysis”, Addison-Wesley, 1999.
[28] M. Petyt, “Introduction to finite element vibration analysis” Cambridge University Press,
New York, USA, 1990.
[29] S. P. Timoshenko, S. Woinowsky-Krieger, “Theory of plates and shells” McGraw-Hill, 1959.
[30] A. W. Leissa, Vibration of Shells, NASA SP-288, 1973.
[31] D. C. Gazis, “Three-dimensional investigation of the propagation waves in hollow circular
cylinders. I. Analytical foundation”, The Journal of the Acoustical Society of America, Vol.
31(5), pp. 568-573, 1959.
[32] D. C. Gazis, “Three-dimensional investigation of the propagation pf waves in hollow circular
cylinders. II. Numerical Results”, The Journal of the Acoustical Society of America, Vol.
31(5), pp. 573-578, 1959.
[33] R. Kumar, R. W. B. Stephens, “Dispersion of flexural waves in circular cylindrical shells”,
Proceedings of the Royal Society of London, pp. 283-297, 1972.
[34] C. R. Fuller, “The effects of wall discontinuities on the propagation of flexural waves in
cylindrical shells”, Journal of Sound and Vibration, Vol. 75, pp. 207-228, 1981.
[35] C. R. Fuller, F. J. Fahy, “Characteristics of wave propagation and energy distributions in
cylindrical elastic shells filled with fluid”, Journal of Sound and Vibration, Vol. 81, pp.
501-518, 1982.
43
[36] R. S. Langley, “Wave motion and energy flow in cylindrical shells”, Journal of Sound and
Vibration, Vol. 169, pp. 29-42, 1994.
[37] R. S. Langley, “The modal density and mode count of thin cylinders and curved panels”,
Journal of Sound and Vibration, Vol. 169, pp. 43-53, 1994.
[38] V. V. Tyutekin, “Helical waves of an elastic cylindrical shell”, Acoustical Physics, Vol.
50(3), pp. 273-277, 2004.
[39] V. V. Tyutekin, “Circumferential and helical normal waves of a cylindrical waveguide:
helical waves in a free space”, Acoustical Physics, Vol. 52(4), pp. 471-476, 2006.
[40] K. Heron, “Curved laminates and sandwich panels within predictive SEA”, Proceedings of
the Second International AutoSEA Users Conference, Troy, Michigan, USA, 2002.
44