B Chouvion

A A Popov

S Mcwilliam

C H J Fox

Vibration Modelling of Complex Waveguide Structures

Keywords: waveguide structures, ray tracing, wave scattering, transmission coefficients, 3D frame structures

In this paper, a general wave propagation approach for the vibration analysis of three-dimensional waveguide structures is presented. The analysis is based on a ray tracing method that considers wave propagation and scattering in the formulation, including the effects of decaying near-field wave components. The waveguide dispersion characteristics of straight and curved beams, which govern the propagation and decay of waves as they travel through the different components, are taken into account for different wave types (flexural, torsional, and axial). Separate elements are coupled using reflection and transmission coefficients, and these are derived for a general connection between several non-planar elements. A matrix formulation is used that offers an automatic and concise method for tackling free and forced vibrations problems of complex ring-beam structural systems. The general methodology is demonstrated on simple applications, and the predictions made using the proposed approach are shown to be in excellent agreement with a conventional finite element analysis, with the advantage of reduced computational costs.

Introduction

Many engineering structures are composed of a number of elements that can be modelled as waveguide components (slender, straight and curved beams). This includes structures like deployable frames in aerospace applications or miniaturised mechanical systems (MEMS) based on the vibrations of sensing elements.

The design and optimisation of such structures could be aided greatly by the availability of an efficient technique for rapidly determining the effects of variations in geometry and dimensions on vibration characteristics such as natural frequencies, mode shapes, and forced response.

The vibrational response of simple structures such as uncoupled beams or rings can be determined without difficulty [1] using approximate solutions such as the Rayleigh-Ritz or modal approaches. For more complex structures, an alternative method is the use of Finite Element (FE) packages; however, they are often computationally expensive if many different FE meshes need to be generated when performing optimisation studies. In contrast, when a method based on exact solutions is used, then not only is the solution valid for all frequencies and no model refinements are required as in the FE method, but the formulation also allows the geometry and size of the structure to be varied easily and efficiently.

To model vibrations of complex three-dimensional (3D) frames using an exact solution, each component is modelled as a single element, clearly decreasing the model complexity compared to FE models, and a general method for combining the individual elements is then needed. Several approaches have been developed which use exact solutions for structural elements; these include the mobility approach [2], the dynamic stiffness method [3], the spectral element method [START_REF] Doyle | Wave propagation in structures[END_REF], and the wave scattering approach [START_REF] Cai | Wave propagation and scattering in structural networks[END_REF].

All of those use a matrix formulation to combine and couple the different elements of the structural network, see e.g. [START_REF] Sack | Matrix Structural Analysis[END_REF] or [START_REF] Livesley | Matrix Methods of Structural Analysis[END_REF]. Langley [START_REF] Langley | Analysis of power flow in beams and frameworks using the direct-dynamic stiffness method[END_REF] used the dynamic stiffness method to analyse power flow in general frames. Banerjee and Williams [START_REF] Banerjee | Coupled bending-torsional dynamic stiffness matrix for Timoshenko beam elements[END_REF] have reviewed various dynamic stiffness methods. The spectral element method uses stiffness matrices established in the frequency domain to describe wave propagation dynamics [START_REF] Doyle | Wave propagation in structures[END_REF]. Combining the advantages of the finite element and spectral methods, a wave-packet procedure [START_REF] Kohno | A finite element procedure for multiscale wave equations with application to plasma waves[END_REF] has been recently presented to solve problems of wave propagation in multi-scale behaviour. The wave scattering approach has been studied extensively by several authors. Mace [START_REF] Mace | Power flow between two continuous one-dimensional subsystems: a wave solution[END_REF][START_REF] Mace | Power flow between two coupled beams[END_REF] determined the power flow between two coupled Euler-Bernoulli beams and the distribution of energy within the system, where each subsystem was one-dimensional; and dynamic analyses of three-dimensional frames and truss-type structures were given by Cai and Lin [START_REF] Cai | Wave propagation and scattering in structural networks[END_REF], and Young and Lin [START_REF] Young | Dynamic response of truss-type structural networks: a wave propagation approach[END_REF]. Beale and Accorsi [START_REF] Beale | Power flow in two-and three-dimensional frame structures[END_REF] developed a matrix method which includes multiple wave modes in order to calculate power flow in 3D frame structures. As the method was intended for high frequency analysis, they used Timoshenko beam theory for the flexural modes. Mei and Mace [START_REF] Mei | Wave reflection and transmission in Timoshenko beams and wave analysis of Timoshenko beam structures[END_REF] also studied free and forced vibration analysis of Timoshenko beams using a wave propagation approach. A similar approach that also considers the scattering and propagation of waves is the ray tracing method [START_REF] Cremer | Structure-Borne Sound[END_REF]. It was employed in [START_REF] Hugin | A physical description of the response of coupled beams[END_REF] with the near-field wave component neglected to predict the response (especially the net transmitted power) of coupled beams subjected to flexural vibrations only.

It was later extended to planar circular curved beam structures including the effects of attenuating components [START_REF] Kang | Free vibration analysis of planar curved beams by wave propagation[END_REF] and to the free vibration analysis of complex planar ring-beam structural systems, including cyclic symmetry [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF]. This matrix-based method offers a compact and systematic methodology that allows complex structures, such as multi-span beams, trusses and aircraft panels with periodic supports to be analysed. The main purpose of this paper is to extend this method to analyse 3D waveguide structures composed of curved and straight beam elements and to provide a general methodology suitable for computer implementation which can be used for free and forced vibration response analyses.

Wave approaches have proved to be powerful for analysing the energy transmission through structural networks [START_REF] Beale | Power flow in two-and three-dimensional frame structures[END_REF], and this characteristic has been used to calculate support losses in MEMS structures [START_REF] Chouvion | Vibration transmission and support loss in MEMS sensors[END_REF][START_REF] Chouvion | Comparison of different support loss models for MEMS resonators undergoing in-plane vibration[END_REF]. Wave approaches have also demonstrated computational efficiency for mid-and high-frequency ranges [START_REF] Beale | Power flow in two-and three-dimensional frame structures[END_REF].

The manuscript is organised as follows. Section 2 presents the ray tracing approach to calculate free and forced vibration response of waveguide structures.

Section 3 provides the general matrix method for computer implementation and the calculation of the transmission coefficients for the general case of n connected elements. Section 4 applies the proposed method to illustrative examples and applications.

The ray tracing approach

Displacements defined using a wave approach

Neglecting shear deformation and rotary inertia, Euler-Bernoulli theory is used in this paper. The fundamental governing equations for 3D motion in curved or straight beams are well-known [1], and are reproduced in Appendix A. Structural damping can be introduced into the analysis through the use of a complex modulus of elasticity. In this section, the 3D displacements are defined using a wave approach. This approach forms the bases of the ray tracing method, explained in Section 2.2.

Derived from the governing equations of motion, the 3D displacements of curved and straight beams can be expressed as a sum of waves travelling in opposite direction [START_REF] Graff | Wave Motion in Elastic Solids[END_REF]. Waves that define the in-plane motion (in-plane flexural and longitudinal displacements), as well as waves that define the out-of-plane motion (out-of-plane flexural displacement and twisting of the cross-section) are used and characterised independently. It is interesting to note that similarities exist between the in-plane and out-of-plane motion of waveguide structures. For instance in straight beams, the equations governing the twisting of the crosssection and the longitudinal displacement are both second order, whereas the equations governing the in-plane and out-of-plane flexural displacement are both fourth order. In any curved beam element, the in-plane or out-of-plane flexural displacement, denoted v i and v o respectively, and the longitudinal displacement along the centreline or the rotation of the cross-section due to torsion, denoted w i and w o respectively, are given by:

v(s) = X 1 û+ 1 e -ik 1 s + û+ 2 e -ik 2 s +û + 3 e -ik 3 s -X 1 û- 1 e ik 1 (s-L) + û- 2 e ik 2 (s-L) + û- 3 e ik 3 (s-L) , (1a)
w(s) = û+ 1 e -ik 1 s + û+ 2 X 2 e -ik 2 s + û+ 3 X 3 e -ik 3 s +û - 1 e ik 1 (s-L) - û- 2 X 2 e ik 2 (s-L) - û- 3 X 3 e ik 3 (s-L) . (1b)
These general expressions govern the displacements and twisting of the crosssection and are valid for both in-plane and out-of-plane motions. The subscripts 'i' and 'o' relate respectively to in-plane and out-of-plane motion and can be added to all variables in Eqs. (1). û± i (i = 1, 2, 3) is the amplitude of wave type i travelling in the positive or negative s direction; k i is its associated wavenumber. X i represents the ratio of the radial to tangential amplitudes, or the ratio of the axial to twisting amplitudes, and can be calculated from the governing equations of motion, see Appendix A. In Eqs. (1), the waves travelling in the positive s direction (superscript +) originate from the location s = 0, whilst the waves travelling in the negative s direction (superscript -) originate from the location s = L. In applications, it is convenient to choose s = 0 and s = L to be located at each end of the element.

In a straight beam element, the equations governing the displacements are uncoupled and the ratios of wave amplitude become X 1 = 0 and X 2 = X 3 = ∞.

From Eq. (1), the 3D motion of the element (straight or curved Euler-Bernoulli beam) is defined with twelve wave amplitudes, six to characterise the in-plane motion and six to characterise the out-of-plane motion; within each of these, three waves are travelling in each direction. These wave amplitudes are taken as the basic unknowns in the model.

The ray tracing procedure

The ray tracing method has been used to model in-plane free vibration of planar structures in [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF]. Here, it is extended to model 3D (free and forced) vibrations of 3D waveguide structures. The basis of the approach is similar and is briefly presented in this section.

The entire structure is composed of waveguide elements only, see for instance the general ring-beam structural system illustrated in Fig. 1. It is assumed that the motion of each of these elements can be defined as a sum of waves, see Section 2.1. Euler-Bernoulli theory is used here but the method can be extended to Timoshenko theory without fundamental difficulty, by changing the governing equations of motion, the corresponding wavenumbers, and wave amplitudes, see e.g. [START_REF] Beale | Power flow in two-and three-dimensional frame structures[END_REF] or [START_REF] Mei | Wave reflection and transmission in Timoshenko beams and wave analysis of Timoshenko beam structures[END_REF]. The different elements are connected and coupled at discontinuities to form the overall structure. The ray tracing method is a matrix approach, based on exact solutions, similar to the dynamic stiffness matrices approach [3] or the spectral element method [START_REF] Doyle | Wave propagation in structures[END_REF], and the fundamental methodology of this method is to assemble matrices that define the travel and scattering of waves through the structure.

The ray tracing method considers an initial wave amplitude vector a 0 which contains the initial wave amplitude sources in the positive and negative directions for each element. Each element contains wave travelling in both directions, and the starting points of these waves correspond to the element boundaries, see Assuming that the structure is composed of N elements (straight or curved beam), the initial wave amplitude vector has the form:

a 0 = u 1+ initial u 1- initial . . . u N + initial u N - initial T . (2)
Each u j± initial contains the initial wave amplitudes in element j in the positive (superscript +) or negative (superscript -) direction, and has six components (three defining in-plane motion and three defining out-of-plane motion).

Waves travelling from one end of an element to the other end implies amplitude attenuation and phase change, both depending on the wavenumber and the length of the element. These effects are expressed in a diagonal matrix D called the dispersion matrix, whose diagonal elements have the form

D ii = e -ik i L j ,
where k i is the wavenumber associated with wave i and L j is the element length through which it travels.

When a wave impinges on a discontinuity (external boundary or discontinuity between two or more different elements), it scatters in the different elements with a portion transmitted to the other elements and the remainder reflected back. Generally, wave transmission and reflection at a discontinuity introduce wave mode conversion. For instance, a principal flexural wave can give rise to a mixture of flexural, decaying and longitudinal waves. This phenomenon is quantified using complex transmission coefficients that are obtained by applying force equilibrium and displacement continuity at the attachment between the elements, see Section 3.2. With the wave amplitude vector defined by Eq. (2), the wave scattering is expressed by an overall transmission matrix T, where ûi = T ij ûj and T ij is the transmission coefficient from a wave of amplitude ûj to a wave of amplitude ûi .

The wave amplitude vector of the initial trace is a 0 , defined as in Eq. (2).

Consider now the wave amplitude vector a i of the i th trace. The next ray trace is created after the waves have travelled along their respective element and one transmission has taken place. Using the previously defined dispersion matrix D and transmission matrix T, this is expressed as a i+1 = TDa i . The final set of wave amplitudes, vector a, present in the structure, which defines the steady state response, is expressed as the sum of all the previous ray traces, such that

a = a 0 + a 1 + . . . + a ∞ = ∞ i=0 (TD) i a 0 . (3)
The presence of decaying components ensures that the amplitudes of the waves are attenuated from one trace to the next, and |TD| < 1. The terms on the right hand side of Eq. ( 3) form a geometric series with common ratio smaller than one.

Using this fact, Eq. ( 3) can be rewritten as:

a = (I -TD) -1 a 0 , (4)
where I is the [N × N ] identity matrix.

Free and forced vibration analyses

Using the ray tracing method and its fundamental equation (Eq. ( 4)), the free and forced vibration response of any waveguide structure can be derived. Free vibration analysis using the ray tracing method was explained in [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF], whilst the forced response analysis is considered here.

For a forced response analysis, an external force or moment has the effect of injecting waves into the structure, e.g. see [START_REF] Mei | Wave reflection and transmission in Timoshenko beams and wave analysis of Timoshenko beam structures[END_REF]. In the following explanation, the term 'force' represents either an external concentrated force (tensile or shear in-plane and out-of-plane), or an external moment (twisting or bending moment in-plane or out-of-plane), or a combination of these. Only concentrated forces acting on the beams centreline are considered. The excitation point is considered as a particular discontinuity, and separate beam elements are considered on each side of it. The excitation force can also be applied at a junction connecting several elements, such that the external force creates waves in all elements directly attached to the excitation point.

The waves generated by the external forces are taken into account in the a 0 term. For instance, if a force is present between two elements only (the n th and

(n + 1) th ), a 0 becomes:

a 0 = 0 0 1 st element . . . 0 u n- F n th element u (n+1)+ F 0 (n+1) th element . . . 0 0 N th element T . (5)
The terms u n- F and u

(n+1)+ F
represent the waves generated by the force in the negative direction of component n and the positive direction of component (n+1).

These depend on the magnitude of the externally applied forces and can be obtained by considering displacement and slope continuities, and forces equilibrium at the excitation point. Section 3.2 presents a general approach to calculate these waves, depending on the junction characteristics.

Substituting the value of a 0 into Eq. ( 4), the final set of wave amplitudes a can be calculated by solving the resulting equation. These wave amplitudes can be used in Eqs. ( 1) to calculate the displacement response of the structure at the excitation frequency considered.

For free analysis, there are no initial amplitudes in the structure and the term a 0 contains only zeros. In this situation, Eq. ( 4) becomes

(I -TD)a = 0. ( 6
)
This equation is similar to that used in the wave train closure principle [START_REF] Cremer | Structure-Borne Sound[END_REF].

The natural frequencies of the structure can be obtained by solving numerically the following characteristic equation derived from Eq. ( 6)

|I -TD| = 0. ( 7
)
The corresponding mode shapes can be determined by back-substituting the frequency solutions in the T and D matrices, which are both frequency dependent, and calculating the corresponding wave amplitude vector a by solving Eq. ( 6).

These wave amplitudes can be used in Eqs. (1) to determine the mode shapes.

Implementation

The ray tracing method, as presented in Section 2, is a general method to analyse free and forced vibration of waveguide structural systems. The following section presents a particular formalism to define the geometry of complex structural networks.

Definition of structure geometry

In order to study any complex waveguide system, such as the one shown in Fig. 1, the overall geometry of the structure must be specified from knowledge of the position and coupling between the different elements (straight or curved beams). The displacements and wave amplitudes of each element are defined in a local coordinate frame. For curved beam elements, the associated local coordinate frame corresponds to a circumferential coordinate system.

All elements are numbered from 1 to N . Each element j and its associated local coordinate frame (x, y, z) j , is defined relative to another element to which it is connected. Consider the case when the k th element is connected directly to the j th element, the connection between element j and element k is defined with a rotation of the j th coordinate system with respect to the k th coordinate system. This rotation is expressed in terms of the vector direction of the axis of rotation and the associated angle of rotation. With a normalised vector direction

u x u y u z j
of the axis of rotation, expressed in (x, y, z) j , and an angle of rotation α about this axis, then the transformation matrix of rotation M k→j rot is:

M k→j rot =      u 2 x (1 -c) + c u x u y (1 -c) -u z s u x u z (1 -c) + u y s u x u y (1 -c) + u z s u 2 y (1 -c) + c u y u z (1 -c) -u x c u x u z (1 -c) -u y s u y u z (1 -c) + u x c u 2 z (1 -c) + c      , (8)
where c = cos α and s = sin α. This matrix defines the rotation from frame k to frame j, such that:

U j k = M k→j rot U k k , (9)
where U j k contains the displacements (or slopes, forces, moments) of the k th element expressed in the local frame j, and U k k contains the corresponding displacements (or slopes, forces, moments) of the k th element expressed in its local frame k.

By assuming that the local coordinate frame of one particular element corresponds to the global coordinate frame, then all the different elements can be expressed in the global coordinate system by the successive uses of Eq. ( 9). If curved beam elements are present, additional matrices of rotation must be used to simulate the rotation of the local frame along the corresponding element.

To create the geometry of the entire structure, all the straight beam elements are defined with the following: element number, material properties, cross-section properties, length, and rotation characteristics (vector direction and angle) with respect to a particular element on which it is directly connected. In the case of curved beam elements, one must add its radius of curvature and the 'length' corresponds to the circumferential length. It is also important to notice that for curved beam elements, the associated local frame is expressed with respect to the curvilinear position along the element and the local frame actually rotates in the global system. Boundary conditions can also be added to the model by specifying their type and location (element number, and positive or negative end of the element). This particular formalism is however not essential for the proposed wave approach and other types of geometry definition can be used. The approach described above is not limited to simple structures composed of beams connected at various angles. Indeed, several elements (straight or curved beams) can be connected at the same point.

The connection between several elements, expressed using a rotation matrix, gives rise to particular transmission coefficients. These can be calculated using the method explained in the following section.

Transmission coefficients calculation

The application of the ray tracing method relies on having knowledge of the detailed transmission characteristics of waves at discontinuities of the structure.

The reflection and transmission of waves in Euler-Bernoulli beams for different discontinuities have been determined by Mace [START_REF] Mace | Wave reflection and transmission in beams[END_REF], and by Chouvion et al. [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF] for a planar ring/beam transmission. This section presents a general method to derive the transmission coefficients of n beam elements attached at a joint. These beam elements can be straight or curved, not necessary in the same plane, and are defined using the wave approach explained in Section 2.

The transmission coefficients are calculated by considering the continuity and force equilibrium equations at discontinuities taken in isolation from the rest of the structure. Each joint is studied independently. It is possible to introduce energy loss for non-conservative discontinuities in which the total power impinging and leaving the discontinuity is different. However, this is not considered here.

External forces, or moments, acting on the joint can also be included. The reflection coefficients at boundaries can be found without difficulty by using the equations governing the boundary conditions, see e.g. [START_REF] Cremer | Structure-Borne Sound[END_REF][START_REF] Graff | Wave Motion in Elastic Solids[END_REF] for straight beams, or [START_REF] Kang | Free vibration analysis of planar curved beams by wave propagation[END_REF] for curved beams.

The model considers the neutral axes of the elements to be connected at a single point. The dimensions of the joint connecting the beams and the thickness of the beams are assumed to be negligible. Flexibility of the joint might be included, typically by matching a finite element model of the joint to wave motion in the beams. The transmission coefficients for two beams joined at an angle, taking into account the geometry and mass of the joint, have been considered in [START_REF] Mei | Wave reflection and transmission in Timoshenko beams and wave analysis of Timoshenko beam structures[END_REF]. In addition, the transmission between a ring and a beam has been considered in [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF]. A general calculation for the transmission coefficients for n elements (straight or curved beams) connected in a 3D structure is presented here. The method is explained in this section and uses the transformation matrix defined by Eq. ( 8).

The wave incident on a joint along a beam element, which can be straight or curved, can be of any type (predominantly longitudinal, in-plane flexural propagating, etc.). The presence of the joint ensures that part of this wave is reflected back along the same element and the remainder is transmitted into the other elements. This partial reflection is also accompanied by mode conversion, so that the incident wave generates other types of wave in each of the elements. The local transmission matrix T n , of a joint connecting n elements, contains the transmission coefficients from any wave type, in any element, to any resulting wave after transmission or reflection; and is of size [6n × 6n].

Assembly of the equilibrium and continuity expressions at the joint yields a system of 6n simultaneous equations that can be solved to provide values for the required transmission coefficients for each wave type. In the following description, the joint connects beam elements numbered from 1 to n. It is assumed that the joint is situated at the positive end of the element whose frame is the frame of reference (element 1 in this case) and at the negative end of the others elements (2, . . . , n). In each element, the displacements are defined as a sum of waves (see Section 2.1) travelling in the positive and negative directions. External forces and moments acting on the joint are also included.

The displacements and slopes (or twisting) continuities between element 1 and element j (j = 2, . . . , n) give the following equations:

X 1 = M j→1 rot X j , (10a)
Φ 1 = M j→1 rot Φ j , (10b)
where X j and Φ j are vectors containing the displacements and rotations of the component j, respectively, in the x, y and z local directions, taken at the attach-ment point. They are defined as:

X j = v j i v j o w j i T , (11)
Φ j = - ∂v j o ∂s φ j i w j o T , (12)
where φ j i is the rotation of the cross-section along the y axis. For straight beams,

φ j i = ∂v j i ∂s
, whilst for curved beams φ j i =

∂v j i ∂s + w i R [1].
The forces and moments of each element are related as follows:

F 1 = n j=2 M j→1 rot F j + F ext , (13a)
M 1 = n j=2 M j→1 rot M j + M ext , (13b)
where F j and M j contain the internal forces and moments of element j along the

x, y and z directions, respectively. F ext and M ext are the applied external forces and moments, defined in the coordinate system linked to element 1. The force and moment vectors considered are defined by three components:

F = S P T T , (14)
M = M 2 M 1 M T T . ( 15
)
where S, P and T are the in-plane shear force, the out-of-plane shear force and the normal force, respectively; and M 2 , M 1 and M T are the out-of-plane bending moment, in-plane bending moment, and twisting moment, respectively.

Eqs. ( 10) and ( 13) give 6n equations that can be solved simultaneously to find the required transmission coefficients. Consider an incident set of waves a incident , a transmitted/reflected set of waves a created , and a transmission coefficient matrix

T n defined as:

a created = u 1-u 2+ u 3+ . . . u n+ T , (16)
a incident = u 1+ u 2-u 3-. . . u n- T , (17)
and

T n =         B 1+ Π 2→1 . . . Π n→1 Π 1→2 B 2-. . . Π n→2 . . . . . . . . . . . . Π 1→n Π 2→n . . . B n-         , (18)
such that

a created = T n a incident , (19)
where Π i→j contains the transmission coefficients from waves incident from element i to element j. When i = j, the transmission coefficients are reflection coefficients denoted by B i± , where ± indicates the positive or negative end of the component. Eqs. ( 10) and ( 13) can be combined and re-arranged in the form

T n = C -1 H (20)
to calculate the transmission matrix T n , where C and H are [6n × 6n] matrices.

The transmission coefficients for two beams connected at an angle and lying in a plane in their local coordinate system are presented in Appendix B.

To calculate the amplitudes of the waves generated by externally applied forces and moments, Eqs. ( 10) and ( 13) can be combined and re-arranged in the form

u 1- F u 2+ F u 3+ F . . . u n+ F T = C -1      0 [6(n-1)×1] F ext M ext      . ( 21
)
Appendix C considers the particular case when forces and moments are applied to a straight beam.

The transmission coefficients for each discontinuity can be found using Eq. ( 20).

To apply the ray tracing method, the local transmission matrices are assembled in the global matrix T and Eq. ( 4) is solved as described in Section 2.3. Some examples illustrating the application of the ray tracing method are presented in Section 4.

Examples

Transverse force at the end of a cantilever beam

To illustrate the forced response analysis using the ray tracing method, presented in Section 2, the simple case of a cantilever beam of length L subjected to an externally applied transverse force is presented. The system studied is shown in Fig. 2. Only in-plane flexural vibrations are considered (in the (x, z) plane), and the number of wave amplitude unknowns needed to define the displacements is four per beam element (flexural propagating and decaying near-field waves travelling in each direction).

To consider the general case when the externally applied force acts anywhere along the length of the beam, the beam is split into two elements. Element 1 has one end at clamp and the other at the excitation point, and has length L 1 .

Element 2 has one end at the excitation point and the other at the free end, and

has length L 2 = L -L 1 .
The wave amplitude vector is given by:

a = û1+ 2 û1+ 3 û1- 2 û1- 3 û2+ 2 û2+ 3 û2- 2 û2- 3 T , (22)
and the corresponding initial wave amplitude vector created by the applied force (see Appendix C) is:

a 0 = S ext 4EI y k 3 Fy 0 0 -i -1 -i -1 0 0 T , (23)
where S ext is the applied force, E is the material Young's modulus, I y the second moment of area of the cross-section, and k Fy is the wavenumber for flexural waves.

The transmission matrix T [8×8] and the dispersion matrix D [8×8] are defined as:

T =         0 T clamped 0 0 0 0 0 I I 0 0 0 0 0 T free 0         , (24)
D =         ∆ L 1 0 0 0 0 ∆ L 1 0 0 0 0 ∆ L 2 0 0 0 0 ∆ L 2         , (25)
where

T clamped =   -i -1 -i -1 + i i   and T free =   -i 1 + i 1 -i i 
 are the reflection coefficient matrices at the clamped and free boundaries respectively [START_REF] Graff | Wave Motion in Elastic Solids[END_REF]; and ∆ L j (j = 1, 2) is defined for in-plane flexural vibration as

∆ L j =   e -ik Fy L j 0 0 e -k Fy L j   .
Solving Eq. ( 4) with the above appropriate expressions for a, a 0 , T and D, and with L 1 = L (force applied at the tip of cantilever); and substituting the obtained wave amplitudes into Eq. ( 1), the following exact expression for the transverse displacement w i of the beam can be derived:

w i (z) = S ext 2EI y k 3 Fy Φ e ik Fy z -Θ e -ik Fy z + Γ e k Fy z -Ψ e -k Fy z ∆ , (26)
where:

Φ = f 1 + i + 2if g + (i -1)g 2 , Θ = (i -1)f + (i + 1)f g 2 + 2ig, Γ = g 1 + i -(i -1)f 2 + 2f g , Ψ = (1 -i)g + (1 + i)f 2 g + 2f, ∆ = 1 + f 2 + 4f g + f 2 g 2 + g 2 ,
f = e -ik Fy L , and g = e -k Fy L .

It can be seen from Eq. ( 26) that the displacement of the beam (flexural modes only) is infinite when ∆ = 0. This condition indicates the presence of natural frequencies, and is the standard expression for the natural frequencies of a clamped-free beam, see e.g. [1].

Free vibration of a ring-based resonator

The MEMS ring-based resonator shown in Fig. 3 has been studied by several authors [START_REF] Hopkin | Performance and design of a silicon micromachined gyro[END_REF][START_REF] Fell | A second generation silicon ring gyroscope[END_REF] and exhibits cyclic symmetry. In the study performed in [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF],

based on the work by Thomas [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF], it was shown that the response of each (cyclically symmetric) sector can be linked to its neighbours by a particular change of phase, such that the in-plane displacement on one side is equal to e -i 2π N κ times the in-plane displacement on the other side; where N is the number of identical sectors that form the overall structure, and κ is the cyclic mode number. In the ray tracing method this is achieved by incorporating the change of phase as a transmission coefficient. The same cyclic symmetry simplification is used in the present study, and the entire structure (Fig. 3(b)) is modelled using only one sector consisting of three connected beams and a portion of the ring. As the structure lies in a plane (the (x, z) plane), the direction vectors used to define the connection between the different beam elements, see Section 3, are all equal to 0 1 0 . To study both in-plane and out-of-plane motions simultaneously, twelve wave amplitude unknowns per elements are needed, which make a total of sixty unknowns. It is assumed that the cross-sections of the ring and legs are uniform, and that the material (silicon) is also uniform. This ensures that the in-plane and out-of-plane modes of the resonator are uncoupled and can be analysed separately.

The material properties used are: mass density ρ = 2329 kg/m 3 , E = 170 GPa, and Poisson's ratio ν = 0.28. The structure dimensions, same as those given in [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF], are as follows. The clamped beam, length 1.047 mm, is connected at an angle 94.5 • to a second beam, length 2.582 mm, itself connected at an angle 54.5 • to a third beam, length 1.064 mm. This third beam is attached to the ring with radius 2.94 mm. The in-plane width of the entire leg is 0.061 mm, whilst the inplane width of the ring is 0.121 mm. The out-of-plane (axial) width of the entire structure is 0.1 mm.

Five different analyses (κ = 0, 1, . . . , 4, see [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF]) are performed to calculate the first few natural frequencies for each value of κ. The calculated frequencies are compared with FE predictions using Euler-Bernoulli and Timoshenko two-dimensional (2D) beam elements, as well as a 3D model. The latter model gives a more accurate representation of the resonator by using curved segments to connect the two straight beam elements in the legs and modelling the fillet radii used in practice to avoid stress concentrations at the joints. This is achieved by meshing the structure with twenty-noded 3D brick elements. Results for the in-plane and out-of-plane natural frequencies are shown in Table 1. Excellent agreement (up to five significant digits), for low and higher order modes, is obtained with the ray tracing method and the FE analysis using Euler-Bernoulli elements. The theory behind these methods is identical, where only the element centrelines are modelled, and shear deformation and rotary inertia are neglected.

With Timoshenko beam elements, the differences are still small (less than 0.5%) and generally increase as the mode number increases, due to shear deformation effects for higher order modes. The differences between results from the ray tracing method and the 3D finite element analysis are slightly larger. This may be explained by the presence of 'rounded' corners, and also as the number of elements used (twenty along the length of each component, two along their width, and two along their thickness) is probably insufficient to accurately model flexure and torsion of the curved beams. A refined mesh gives results for the natural frequencies slightly closer to the ones obtained using the ray tracing method.

The lowest modes with κ = 2 and κ = 3 (called respectively 2θ and 3θ modes) are the modes used in practical applications of the device [START_REF] Fell | A second generation silicon ring gyroscope[END_REF]. Wong [START_REF] Wong | Thermoelastic damping in MEMS ring resonators[END_REF] performed experiments on such a resonator and measured natural frequencies for the 2θ and 3θ modes with the use of a laser vibrometer. These measured values are compared with those obtained using the ray tracing method in Table 2. Wong [START_REF] Wong | Thermoelastic damping in MEMS ring resonators[END_REF] developed an analytical model for the vibrations of the ring-based resonator. The approach used the static deflection of the leg to estimate the strain and kinetic energies, and was based on the assumption that the displacement of each leg at the point of attachment to the ring is identical to that of the free ring. It was also assumed that the presence of the legs did not significantly affect the mode shape of the ring. The results shown in Table 2 indicate that this model reasonably accurately predicted the in-plane natural frequencies compared to FE results. However, the agreement for out-of-plane natural frequencies is much less satisfactory. In contrast, the ray tracing method is an exact solution, based only on the assumption that ring and leg can be modelled as a one-dimensional waveguide, and the method shows much better agreement with the FE model and experiments.

The mode shapes obtained using the ray tracing method for the lowest natural frequency with κ = 0, 1, . . . , 4 are shown in Fig. 4, and are compared to those found using the FE model using brick elements. Two orthogonal modes exist for κ = 1, 2, 3 [START_REF] Chouvion | In-plane free vibration analysis of combined ring-beam structural systems by wave propagation[END_REF] but only one of these is shown in Fig. 4 for simplicity. The first mode shape with κ = 0 corresponds to rigid body mode of a free ring rotating about its centre for in-plane motion, and a translation of the ring along the axis of symmetry for out-of-plane motion. The first mode shape with κ = 1 corresponds to in-plane translation or out-of-plane rotation (tilting) of the ring. The mode shapes obtained with both methods generally show a good agreement. However, small differences still occur because the wave approach neglects shear deformation and rotary inertia effects. Also, in the ray tracing method, it is assumed that the beams making up the legs are connected at abrupt angle and that each leg is connected at a single point to the ring; this is not the case for the FE model (with brick elements) and the device measured by Wong. This is expected to have little influence on the natural frequencies, but it may have a stronger influence on the modes shapes, especially the leg deformation.

Forced response analysis of a complex 3D structure

The aim of this example is to illustrate the analysis for a complex 3D structure that contains straight and curved beam elements, and to show how the matrices are assembled. The system studied is illustrated in Fig. As explained in Section 3, the connection between the different elements j and (j + 1) are defined using the angle of rotation (and associated vector direction), from the local frame j to the local frame (j + 1). Before defining the different rotations between the elements, it is firstly important to notice that the equations used (and especially the sign convention) that govern the motion of the curved beam elements, are only valid when the direction x is radial inwards, which means that particular attention needs to be given when defining the local frame linked with curved beam elements. With the local coordinate frames, defined as in Fig. 5, the angle of rotation α i→j and associated vector direction v i→j between element i and element j are expressed as:

α 1→2 = α 4→5 = α 6→7 = 2π/3, (27a)
v 1→2 = v 4→5 = v 6→7 = √ 3/3 [-1, 1, 1], (27b)
α 1→3 = α 3→4 = α 4→6 = α 6→1 = 0. ( 27c
)
For indication, v i→j is expressed in the local frame of the i th element.

Using a wave amplitude vector a defined as:

a = u 1+ u 1-. . . u 7+ u 7- T , (28)
the overall transmission matrix T can be assemble as follows. The size of T is

[84×84]. However, all its components can be grouped in blocks as [6×6] matrices, defining the waves travelling in a single element in one direction. T contains null block-matrices except in the following places:

T(1, 2) = B 1-, T(1, 11) = Π 6→1 , T(1, 14) = Π 7→1 , (29a)
T(2, 1) = B 1+ , T(2, 4) = Π 2→1 , T(2, 6) = Π 3→1 , (29b)
T(3, 1) = Π 1→2 , T(3, 4) = B 2-, T(3, 6) = Π 3→2 , (29c)
T(4, 3) = T clamped , (29d)
T(5, 1) = Π 1→3 , T(5, 4) = Π 2→3 , T(5, 6) = B 3-, (29e)
T(6, 8) = I, (29f)
T(7, 5) = I, (29g)
T(8, 7) = B 4+ , T(8, 10) = Π 5→4 , T(8, 12) = Π 6→4 , (29h)
T(9, 7) = Π 4→5 , T(9, 10) = B 5-, T(9, 12) = Π 6→5 , (29i)

T(10, 9) = T clamped , (29j)
T(11, 7) = Π 4→6 , T(11, 10) = Π 5→6 , T(11, 12) = B 6-, (29k)
T(12, 2) = Π 1→6 , T(12, 11) = B 6+ , T(12, 14) = Π 7→6 , (29l)
T(13, 2) = Π 1→7 , T(13, 11) = Π 6→7 , T(13, 14) = B 7-,

T(14, 13) = T clamped .

In Eqs. (29), the matrix T clamped is the [6 × 6] matrices defining the reflection coefficients at a clamped boundary [START_REF] Chouvion | Vibration transmission and support loss in MEMS sensors[END_REF]. If all of the straight beam elements have the same material properties and cross-sections, the structure exhibits cyclic symmetry. In this situation, similar simplification to the one performed in Section 4.2 can be applied to reduce the size of the matrices involved in Eq. ( 4) to [36 × 36].

From the symmetry of the system, some of the transmission coefficients are identical; for instance,

Π 1→2 = Π 4→5 = Π 6→7 , or B 1+ = B 4+ = B 6+ , etc.
To illustrate how the assembly of the overall T matrix is performed and how Eqs. (29) can be obtained automatically with an algorithm, consider Eq. (29b).

Eq. (29b) defines the second line of T, and from Eq. (28), this second line shows how the set of wave amplitudes u 1-is created. These amplitudes are simply created by: a reflection of the set u 1+ at the end of element 1 (u 1+ is in the first position of a -see Eq. ( 28), so these reflection coefficients are placed in T(2, 1)), and transmission from incident sets of waves u 2-and u 3-(fourth and sixth position in a respectively, and so are taken into account in T(2, 4) and T(2, 6), respectively).

The overall diagonal dispersion matrix D is simply defined as:

D =            ∆ L 1 0 . . . 0 0 0 ∆ L 1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 . . . ∆ L 7 0 0 0 . . . 0 ∆ L 7            , (30)
where element j has length L j . For the curved beam elements, L 1 = L 6 = 2πR/3 and L 3 = L 4 = πR/3, where R is the radius of curvature of the ring. Each ∆ L j is a diagonal matrix of size [6 × 6] with diagonal elements:

∆ L j (i, i) = e -ik i i L j , ∆ L j (i + 3, i + 3) = e -ik o i L j , (31)
with i = 1, 2, 3.

With an external force applied between elements 3 and 4, waves are created corresponding to u 3-and u 4+ (sixth and seventh entry in a respectively). On this basis, a 0 has non-zero elements in rows 6 and 7.

With the geometry and matrices defined in Eqs. ( 27)-(30), the ray tracing method is employed to calculate the displacement in three directions (u x , u y and u z , in the radial, axial and tangential directions, respectively) at the point A induced by a concentrated force, acting between elements 3 and 4, as shown in Fig. 5. To excite all modes, the external force is such that: S ext = P ext = -T ext = 20 N, where S ext , P ext and T ext are respectively along the x, y, and z directions of the local frame 3. This response is compared to the one obtained with a FE model in Fig. 6. As can be seen in the figure, there is a very good correlation between results from the ray tracing method and the FE model.

Conclusions

The method presented predicts the free and forced vibration response of a structure composed of straight and curved beam elements. The method is based on a ray tracing approach, which considers waves travelling and scattering through the different structural members. The elements are coupled together with the use of transmission coefficients and these were derived for a general connection between several components. The outcome is a simple matrix formulation which permits a relatively straightforward implementation, and with significant computational efficiency. The method provides an accurate calculation of the response based on exact solutions, but at the same time is more general (travelling waves) that the dynamic stiffness method (standing waves). This has not been exploited in the current paper but the potential has been outlined.

Euler-Bernoulli beam theory was used in this paper and the extension to Timoshenko beams would be an useful contribution for studying structure containing thicker beams and/or for high frequency analysis. The method would be similar to the general method presented here, and only the basic governing equations of motion are expected to change. For a practical study of aerospace frame structures that generally exhibit periodically symmetric properties, the incorporation of periodic symmetry simplification in the formulation will be assessed in future publications.

defining in-plane motion can be expressed as [2]:

-EI y ∂ 3 ∂s 3 w i R + ∂v i ∂s + EA R ∂w i ∂s - v i R = ρA ∂ 2 v i ∂t 2 , (A.1a) EI y R ∂ 2 ∂s 2 w i R + ∂v i ∂s + EA ∂ ∂s ∂w i ∂s - v i R = ρA ∂ 2 w i ∂t 2 ; (A.1b)
whereas the equations defining out-of-plane motion can be expressed as [2]:

GJ R ∂ 2 ∂s 2 w o + v o R + EI x ∂ 2 ∂s 2 w o R - ∂ 2 v o ∂s 2 = ρA ∂ 2 v o ∂t 2 , (A.2a) GJ ∂ 2 ∂s 2 w o + v o R - EI x R w o R - ∂ 2 v o ∂s 2 = ρI zz ∂ 2 w o ∂t 2 . (A.2b)
The moments M 2 , M 1 and M T , in the x, y and z directions respectively, are defined by:

M 2 = EI x w o R - ∂ 2 v o ∂s 2 , (A.3) M 1 = EI y ∂ ∂s ∂v i ∂s + w i R , (A.4)
M T = GJ ∂ ∂s w o + v o R . (A.5)
The forces S, P , and T , in the x, y and z directions respectively, are defined by:

S = - ∂M 1 ∂s , (A.6) P = ∂M 2 ∂s + M T R , (A.7) T = EA ∂w i ∂s - v i R . (A.8)
For a curved beam element, the radial and tangential displacements are coupled, so are the axial displacement and twisting of the cross-section. The ratios defining this coupling can be found from the governing equations of motions (A.1)

for in-plane motion, and (A.2) for out-of-plane motion. These can be calculated as [3]:

X i i = iERk i i I y k 2 i i + A ρR 2 Aω 2 -EI y R 2 k 4 i i + EA , (A.9)
for in-plane motion, and

X o i = k o i R 2 GJ + EI x ρR 2 Aω 2 -EI x R 2 k 4 o i + GJk 2 o i , (A.10)
for out-of-plane motion. In these expressions ω is the circular frequency, k i i and coefficient matrix defined as in Eqs. ( 16)-( 18), Eq. ( 20) gives: Table 2: In-plane (IP) and out-of-plane (OOP) natural frequencies (kHz) for the ring-based resonator obtained with the ray tracing method (RT), a FE model using brick elements, experiments, and a simplified analytical model [START_REF] Wong | Thermoelastic damping in MEMS ring resonators[END_REF].

k o i (i = 1, 2 ,
C =                                  0
                                 , (B.1)
H =                                  0
                                 , (B.2)

Figures and tables

Fig. 1 .

 1 Fig. 1.

5 .

 5 It consists of a ring supported by three clamped straight beams (columns) equally spaced. The straight beam elements are numbered 2, 5 and 7, and the curved beam elements (ring arcs) are numbered 1, 3, 4 and 6. The connection between the different elements follows what is shown in Fig. 5. An external force is applied at the connection point between elements 3 and 4, and explains why that portion of the ring was split into two different elements. Structural damping is assumed to be present in the structure and is taken into account using a complex modulus of elasticity.

  3) are respectively the wavenumbers for in-plane waves and outof-plane waves in the element, obtained by solving the governing equations of motion.For a straight beam element, Eqs. (A.1)-(A.10) can be simplified by letting R → ∞, and the well-known equations of motion for a straight beam are obtained.

Figure

  Figure 1: Complex waveguide structure composed of curved and straight beam elements. Each element has an associated coordinate frame, and contains waves travelling in both directions.

1 :

 1 Figure 1: Complex waveguide structure composed of curved and straight beam elements. Each element has an associated coordinate frame, and contains waves travelling in both directions.

Figure 2 :

 2 Figure 2: Cantilever beam excited by an external transverse force.

Figure 3 :

 3 Figure 3: (a) Photograph of a ring-based resonator on its silicon wafer (courtesy Atlantic Inertial Systems); (b) Schematic representation.

Figure 4 :

 4 Figure 4: In-plane and out-of-plane mode shapes of the ring-based resonator for the lowest natural frequencies of each value of κ, obtained with the ray tracing method (first and third rows) and with a 3D FE model using Abaqus TM (second and fourth rows).

Figure 6 :

 6 Figure 6: Frequency response function at point A calculated with the ray tracing method ('• • • •' ux, '-• -•' uy and '--' uz) and a FE model ('××' ux, '++' uy and 'o o' uz) with ten Euler-Bernoulli beam elements per component (straight or curved). The material used is steel, with ρ = 7860 kg/m 3 , E = 207(1 + iη) GPa, η (loss factor)= 0.02, ν = 0.3. Square cross-sections are used with a width of 10 cm for the straight beams, and 5 cm for the ring arcs. The lengths of the different elements are defined with L2 = L5 = L7 = 1 m and R = 0.8 m.

Figure A. 1 :

 1 Figure A.1: 3D motion of a curved beam: notation and sign convention.

Table 1 :

 1 In-plane (IP) and out-of-plane (OOP) natural frequencies (kHz) for the ring-based resonator calculated with the ray tracing method and different FE models: 2D Euler-Bernoulli and Timoshenko (T'nko) beams (shear correction factor 0.85) elements, and 3D brick elements.

			Ray Tracing	FE analyses 2D Euler 2D T'nko 3D brick
	κ = 0	IP OOP	18.53 40.89 4.820 30.12	18.53 40.89 4.820 30.12	18.48 40.83 4.817 30.07	19.31 42.94 5.025 31.49
	κ = 1	IP OOP	11.96 38.68 5.431 37.29	11.96 38.68 5.431 37.29	11.93 38.61 5.429 37.21	12.61 40.85 5.685 39.05
	κ = 2	IP OOP	14.27 38.08 10.83 46.38	14.27 38.08 10.83 46.38	14.25 38.01 10.82 46.26	14.58 40.30 11.21 48.00
	κ = 3	IP OOP	32.26 41.23 22.91 52.59	32.26 41.23 22.91 52.59	32.21 41.13 22.89 52.43	32.79 43.70 23.65 53.52
	κ = 4	IP OOP	35.88 55.28 30.01 58.86	35.88 55.28 30.01 58.86	35.83 55.06 29.95 58.66	37.25 58.18 31.18 59.70

Acknowledgements

The authors gratefully acknowledge the support for this work provided by Atlantic Inertial Systems.

Appendix A: Governing equations of motion, and forces and moment definitions Consider a curved beam in the (x, z) plane having constant radius of curvature, vibrating in a 3D motion, as shown in Fig. A.1. The curved beam is characterised by the radius of curvature R; the cross-sectional area A; the second moments of area of the cross-section I x and I y in the x and y directions, respectively; the polar moment of inertia I zz ; the Young's modulus E; the mass density ρ; and, the torsional rigidity GJ (where G is the shear modulus and J is often called the torsion constant [1]). I zz and J are identical for circular cross-sections, but for other shapes, J must be determined by other means [1]. In Fig. A.1, S is the in-plane shear force (radial direction), P is the out-of-plane shear force (axial direction), T it the normal force (tangential direction), M 2 is the out-of-plane bending moment (for rotation about x), M 1 is the in-plane bending moment (for rotation about y), and M T is the twisting moment (for rotation about z). v i and v o are the displacements along the x and y axes, respectively; whereas w i and w o are the tangential displacement (along the z-axis) and the angular rotation of the cross-section due to torsion, respectively. The circumferential coordinate along the centreline is denoted by s.

Neglecting the effect of rotary inertia and shear deformation, the equations Appendix B: Special case of the transmission between two beams in the same plane Consider two beams having the same material properties, connected at an angle α, such that the system lies in the same plane (x, z). For this particular case, and with the vector of incident waves, vector of created waves, and transmission where γ = ω ρI y EA , δ = i GJI zz EAI x , s = sin α, c = cos α, and k Fx is the out-ofplane flexural wavenumber.

Substituting these expressions for C and H into Eq. ( 20) gives the required transmission coefficient matrix. The matrix C can also be used in Eq. ( 21) to calculate the waves generated by external forces and moments applied at the junction. In the case where α = 0, solving Eq. ( 21) gives the waves generated expressed in Appendix C.

Appendix C: Waves generated by concentrated forces and moments in a straight beam

Consider the case of a straight beam subjected to external concentrated forces and moments such that the beam is split into two elements, one on each side of the point of application. For this special case, the transmission between element 1 (left side) and 2 (right side) is modelled using Eqs. ( 10)-( 21) with α = 0 and n = 2. These equations can be easily solved analytically, giving a transmission matrix T 2 equal to the identity matrix as the waves impinging on one side of the disturbance travel without interference to the other side.

Concerning the waves generated by the applied forces and moments, the following notation that considers separately in-plane and out-of-plane motion is used

where the wave amplitude vector components are taken in the particular order: mainly longitudinal (or twisting -for the out-of-plane case), mainly flexural propagating, and mainly flexural decaying, such that gives for a straight beam: