Analysis of vibration and wave propagation in cylindrical grid-like structures

The wavepropagation in and the vibration of cylindrical grid structures are analyzed. The grids are composed of a sequence of identical elementary cells repeating along the axial and the circumferential direction to form a two-dimensional periodic structure. Two-dimensional periodic structures are characterized by wavepropagation patterns that are strongly frequency dependent and highly directional. Their wavepropagation characteristics are determined through the analysis of the dynamic properties of the unit cell. Each cell here is modelled as an assembly of curved beam elements, formulated according to a mixed interpolation method. The combined application of this Finite Element formulation and the theory of two-dimensional periodic structures is used to generate the phase constant surfaces, which define, for the considered cell lay-out, the directions of wavepropagation at assigned frequencies. In particular, the directions and frequencies corresponding to wavea tenuation are evaluated for cells of different size and geometry, in order to identify topologies with attractive waveattenuation and vibration confinement characteristics. The predictions from the analysis of the phase constant surfaces are verified by estimating the forced harmonic response of complete cylindrical grids, obtained through the assembly of the unit cells. The considered analysis provides invaluable guidelines for the investigation of the dynamic properties and for the design of grid stiffened cylindrical shells with unique vibration confinement characteristics.


Introduction
Extensive efforts have been exerted in the investigation of wave propagation in periodic structures, which consist of identical elements, or cells, identically connected to each other.Much of the pioneering work on periodic structures was initiated in the middle of the last century by researchers working in the area of solid state physics.A comprehensive description of the initial studies and main results can be found in the 1946 book by Brillouin [1].Mead and his coworkers [2] expanded on Brillouin's work and demonstrated the effects of periodically located impedance mismatch zones in engineering structures such as periodically supported beams, stiffened plates and shells, and grids.An impedance mismatch is generated by discontinuities in geometry and/or material within the structure.From a dynamic point of view, impedance mismatch zones partially reflect incoming elastic waves, thus causing wave interference phenomena.The interference may be constructive or destructive depending on the relative phase between incident and reflected waves.When impedance mismatch zones are periodically located along the structure, destructive interference phenomena occur over specified frequency bands, within which the propagation of waves is impeded.These frequency bands are typically called "stop-bands".This work analyzes cylindrical two-dimensional (2D) periodic structures.In 2D periodic structures, the location and the extension of the stop bands depend on the considered direction of wave propagation on the structure surface [3].This unique directional behavior complements the stop/pass band pattern and makes the application of 2D periodic structures as directional mechanical filters extremely attractive.The analysis of 2D periodic structures is of high practical relevance, since stiffened plates, panels and beam grillages fall in this category.Two-dimensional periodic structures of the kind studied in this paper include orthogonally stiffened cylindrical shells, periodic composite shell structures or grid stiffened composite panels and cylinders.In here, the wave propagation characteristics of grids of various topologies are investigated.The goal of the study is to identify cell lay-outs with attractive wave attenuation characteristics and provide guidelines for the design of grid-stiffened cylindrical components with directional characteristics at desired frequencies.This work extends the previous investigations by Mead and Bardell [6], Langley [3,4], and Ruzzene et al. [5].The considered grids are obtained by the assembly of curved beams arranged on the structure's surface according to rectangular and hexagonal architectures.Hexagonal architectures have demonstrated particularly interesting characteristics in planar 2D structures, as demonstrated by Ruzzene et al. in [5].As in [5], the analysis is carried out by combining the Finite Element (FE) method with the theory of 2D periodic structures as expressed by Bloch theorem [1].The dynamic behavior of the unit cells here is predicted by modelling each beam using curved beam elements, formulated according to the mixed interpolation method [11].The FE model for the unit cell is then used to obtain the phase constant surfaces for the considered cell and to identify directions and frequencies of wave propagation.The predictions from the phase constant surfaces are then verified by evaluating the forced harmonic response of the complete cylindrical structure.The computational cost associated with this operation is significantly reduced by applying the reduction method for rotationally periodic structures proposed by Thomas [7] and Petyt [8].
The paper is organized in 6 sections and one appendix.In Section 1 an introduction is given, while Section 2 summarizes the fundamentals of wave propagation in 2D periodic structures and briefly describes Bloch Theorem and its implications.Section 3 presents the considered curved beam element formulation and presents some examples that validate the element implementation.The methodology used to compute the structure's harmonic response is described in Section 4, while Section 5 presents the results obtained for the considered topologies in terms of phase constant surfaces and harmonic response.Finally Section 6 summarizes the main results of this study and gives some recommendation for future investigations.

Wave propagation in 2D periodic structures
Bloch's theorem [1] is the tool traditionally used to study the wave motion through a periodic system.The theorem states that the proportionate change in wave amplitude occurring from cell to cell does not depend on the cell location within the periodic system.Accordingly, wave motion can be described as follows: In here and in the following, vectors are denoted with underlined greek or latin letters, while matrices are indicated by bold capital latin letters.In Eq. ( 1), x and y define the two dimensions of the structure, w(x, y, n x , n y ) is the displacement of a point x, y belonging to the cell at location n x , n y , while g(x, y) describes the motion of a single cell.Also, µ x and µ y are the propagation constants in the x and y direction.The propagation constants are complex numbers µ k = δ k + iε k (k = x, y), whose real and imaginary parts are denoted respectively as 'attenuation' and 'phase' constants [3].The propagation constants control the nature of elastic wave propagation in a 2D periodic structure.If the propagation constants are purely imaginary, waves are free to propagate, however if a real part exists, attenuation of the wave's amplitude will be observed as it propagates from one cell to the next.In the analysis of wave propagation in 2D periodic structures, the attenuation constants δ x and δ y are typically set to 0, while the phase constants ε x and ε y are varied from −π to π. Solving with respect to frequency for all the combinations of ε x , ε y in the considered range yields a series of functions ω = f (ε x , ε y ).These functions are known as 'phase constant surfaces' and identify the frequencies of free wave motion in the structure for any given pair of ε x and ε y .The phase constant surfaces are generally represented as iso-frequency contour plots in the ε x , ε y plane.It can be shown [3] that the wave's group velocity at a given frequency lies in the direction of steepest ascent on the phase constant surfaces, which can be found as the normal to the corresponding iso-frequency contour line in the ε x , ε y plane.This property not only allows the determination of wave propagation directions in the structure, but also identifies regions within the structure where waves do not propagate.This unique characteristic depends on the geometry of the unit cell and on the frequency of wave propagation.The analysis of phase constant surfaces guides the design of 2D periodic structures where waves at certain frequencies do not propagate in specified directions.
In here, we analyze the characteristics of wave propagation in cylindrical 2D periodic grids, which are obtained by the assembly of curved beams arranged on the structure surface according to rectangular and hexagonal architectures (Figs 1 and 2).The behavior of the unit cell is described through a Finite Element model developed using curved beam elements.Details of the element formulation are given in the next section and in the appendix.The wave propagation analysis is then performed by applying Bloch's theorem to the resulting discretized model of the unit cell, according to the procedure described in Section 4.

Finite element modelling of a unit cell: Curved beam element
The investigation of the considered class of cylindrical grid structures requires the implementation of a general curved beam element.This section discusses the considered formulation for such an element, and presents some examples used for its validation.The formulation of a general curved beam element is primarily complicated by two difficulties: the interpolation of the curved geometry to the beam displacements and the phenomenon of membrane locking [9].The first difficulty is handled by the careful selection of the coordinate and displacement functions so that the resulting Jacobian can relate the curvilinear geometry to the global coordinate system, and the natural  coordinates to the curvilinear geometry.Membrane locking occurs because of the inextensibility condition of the membrane strain when the thickness of the curved beam element approaches zero.Numerical evidence shows that if this phenomenon is not handled appropriately, the resulting element is much too stiff to be of any use.Results as bad as only predicting 5% of the actual displacements are possible if membrane locking is not handled in the formulation.However, in the case of very shallow or straight beam elements, this phenomenon is either negligible or non-existent.
There are several approaches to resolve the issue of membrane locking.A simple way to handle the problem is to use reduced or selective integration when computing the membrane strain energy terms [9].This technique successfully removes the spurious constraints that cause the element to lock, but it can also introduce spurious modes that can lead to inaccuracies.Alternatively, Meck [10] has shown that full integration of quintic shape functions yields a displacement based element that does not exhibit membrane locking.While this is indeed helpful, using a quintic shape function can cause computations to become very intensive.
Bathe's [11] approach is to abandon a purely displacement based formulation, and adopt a mixed interpolation method to derive a simple yet efficient formulation for a general curved beam element.The transverse, bending and membrane strains are related to the nodal displacements and rotations by evaluating the displacement-based strains and equating them to the assumed strains at the Gauss integration points.In other words, the mixed interpolated element matrices can be obtained by numerically evaluating the displacement-based element matrices at the Gauss integration points corresponding to the number of nodes in an element.Bathe's formulation only considers three strains: the membrane strain, and the two in-plane shear strains.The two transverse strains are neglected, while the out of plane shearing is not considered.The formulation however can include torsional warping through the addition of an appropriate out-of-plane warping displacement [12].

Kinematic relations and shape functions
The Cartesian coordinates x, y, z of a point in the element are expressed using the natural coordinates r, s, and t (see Fig. 3):, where, x k , y k , z k are the coordinates of nodal point k, a k , b k are the element's height and width at node k, while h k denote the shape function.Also 0 V k tx , 0 V k ty , 0 V k tz and 0 V k sx , 0 V k sy , 0 V k sz are respectively the components of the unit vectors t and s at node k.
The behavior of the curved beam element is described by 3 displacements and 3 rotations at each nodal location.Accordingly, the vector of the nodal degrees of freedom for node k can be expressed as: The displacement of any point within the element can be expressed in terms of the nodal degrees of freedom using the same functions employed to map the coordinates (Eq.( 2)): and The strain and displacement interpolation matrices can be formed using Eqs ( 2) and ( 7).The displacement interpolation relations are expressed in the following matrix form: where u = {u, v, w} T , and where H is the displacement interpolation matrix, which is found by substituting Eqs ( 5) and ( 6) into Eq.( 4).The nodal displacements in the natural coordinates are related to the displacements in cartesian coordinates through the following expression: where J is the Jacobian matrix.Details of the derivation of the strain-displacement interpolation matrix are given in Appendix A. The resulting relation between strains and nodal degrees of freedom can be expressed as: where B is the strain-displacement interpolation matrix.

Element mass and stiffness matrix
The interpolation of displacements and strains is used to express the element strain and kinetic energy in terms of the nodal displacement vector δ (e) .
The strain energy of a curved beam element can be expressed as: where V is the volume of the beam, and C is the constitutive matrix for the beam material (see Appendix B).Imposing the strain-displacement relations defined in Eq. ( 10) yields the following expression for the strain energy: where K (e) is the stiffness matrix for the considered curved beam element, which is given by: Similarly, the kinetic energy for the curved beam can be written as: where ρ is the density of the beam's material, and where () indicates differentiation with respect to time.Imposing the displacement interpolations defined in Eq. ( 8) gives: where M (e) is the mass matrix for the considered beam element, defined as:

Model validation
The static performance of the general curved beam element is first evaluated.The tip rotation of the curved cantilever beam shown in Fig. 4 is computed using a single 3 node element and is compared with analytical predictions.The comparison is carried out for varying h/R ratios in order to evaluate the element's robustness for increasing radii of curvature for the beam.The properties of the cantilevered beam are listed in Table 1 and the results are listed in Table 2, which illustrates the accuracy of the considered curved beam element.
The dynamic performance of the general curved beam element is validated by comparing the FE predictions with the ANSYS results for a strip of 5 rectangular cells (Fig. 5).The geometry and material properties of the strip are given in Table 3.The first 10 natural frequencies (ignoring the rigid body modes) are found for the FE model and compared with the ANSYS results in Table 4.The results show a very good agreement, which is obtained despite using in the FE model a mesh considerably coarser than the one considered in ANSYS.In ANSYS, the model of the strip had 480 2-node BEAM4 elements, with 476 nodes, while the developed FE model only features 80 3-node general beam elements and 156 nodes.

Bloch reduction
The Finite Element formulation presented in the previous section is here used to relate generalized displacements and forces at the interface between each cell and its neighbors [4].For wave motion at frequency ω, the relation between cell's generalized forces and displacements can be expressed as: where K and M are the assembled mass and stiffness matrices for the cell, and δ, F are respectively the vectors containing generalized nodal displacements and forces of the cell (see Fig. 6): where δ I and F I are nodal displacements and forces internal to the cell.According to Bloch's Theorem, the generalized displacements and forces at the cell interfaces are related by: 20) and ( 21) can be rewritten in matrix form as follows: where  while A and B are matrices defined in Appendix B. Substituting Eq. (20) into Eq.( 17), pre-multiplying the resulting equations for A H , with H denoting a complex transpose conjugate, and assuming F I = 0 gives: where, are stiffness and mass matrices reduced according to Bloch's theorem.Equation ( 25) is an eigenvalue problem whose solution for a given pair of ε x and ε θ gives the frequency of wave propagation ω for the assigned pair of propagation constants.Varying ε x and ε θ from −π to π and assuming the attenuation constants to be zero, generates the phase contour plots for a given cell.The phase constant surfaces are then used to evaluate the characteristics of wave propagation for the structure obtained by the assembly of the cells.Examples of phase constant surfaces for the considered cell topologies will be described in the results section.

Harmonic response of cylindrical grids
The predictions from the analysis of the phase constant surfaces can be evaluated by computing the harmonic response of the complete cylindrical grid at specified frequencies.The computational effort required for such an analysis can be significantly reduced by considering the cylindrical grids as examples of rotationally periodic structures.Rotationally periodic structures consist of a number of identical elements or 'strips' connected to form a closed ring [7].Accordingly, their geometry at a given angular position θ is identical to the geometry at the angle (θ + 2πn/N ), where N is the number of periodic components in the complete structure and n is an integer varying between 0 and N .The analysis of this class of structural systems is simplified by taking advantage of their rotational periodicity, which allows the analysis of the complete structure to be performed by considering a single strip with imposed periodic boundary conditions.

Equivalent nodal load
According to [8], the harmonic response due to a distribution of harmonic forces can be obtained through a series of calculation performed on a single strip.The equivalent nodal force on the s th strip due the applied external loads can be expressed as: where Q p defines the magnitude of the equivalent nodal loads on the p th strip.From Eq. (28), the force on the first component is simply: and the forces in subsequent components are identical except for a phase difference = 2πp/N .Expressing Eq. ( 28) in matrix form gives: where Φ is a N × N Hermitian matrix [7] whose (p, q) element is e (−i2πp(q−1)/N ) .

Reduced equation of motion
The response of the structure to each of the harmonic terms in Eq. ( 28) can be obtained from separate analyses.Each component induces a deflection of the structure with the same frequency and phase variation, which can therefore be expressed as: where δ I and δ L denote vectors containing the intermediate nodal degrees of freedom and those on the left boundary of the considered strip.Matrix W is defined in terms of the phase difference ε associated to the component of the considered set of loads.Matrix W is given by: where R is a rotation matrix that transforms the displacements of the nodes to the right of the considered strip from the global coordinate system into the local coordinates of the left nodes.Such a transformation is required to enforce the compatibility of displacements after imposing the phase difference ε and therefore to satisfy the conditions for rotational periodicity.
The force acting on a single substructure can be considered as composed of a combination of external loads and boundary forces.Accordingly, the total force acting on the s th strip F (s) , can be written as: where the superscripts e and b respectively denote external loads and boundary forces.Using Eq. ( 33) and following a reduction analysis similar to the one presented in Section 4.1, one obtains: The equation of motion of the s th strip can be generally expressed as: where K (s) , M (s) are the respectively the stiffness and mass matrix for the considered strip, while F (s) is the vector of the applied loads.Substituting Eq. (31) in Eq. ( 35) and pre-multiplying by W H gives: where Equation ( 36) is solved N times corresponding to the components in Eq. (28).The complete set of displacements can be then obtained from Eq. (31).For the considered N values of the total displacements are obtained from: where the s th column of ∆ describes the displacement of the s th substructure.

Considered configurations
Both rectangular and hexagonal cell topologies are considered.The performance of cylindrical surfaces obtained through the assembly of each cell configuration is analyzed through the phase constant surfaces and the model reduction technique based on the rotational symmetry of the structure.For both cell topologies, the analysis is performed on a cylinder of radius 1 m and length 10 m.All the analyzed grids are composed of 16 cells around the circumference and 10 cells along the axial direction.The overall dimensions of the cylindrical structure, such as the radius R and the length L, and the material of the beams are also kept constant and equal to the values listed in Table 5.
The rectangular cells (see Fig. 1(a)) are obtained by the assembly of two straight beams and two curved beams of square cross sections.The cross section of the 2 curved beams and the overall dimensions of the cell are kept constant, while the cross section of the straight beam (thick lines in Fig. 1(a)) is varied in order to evaluate its influence on the wave propagation characteristics of the assembled structure.The cross section for the curved beams measures a = 0.01 m and b = 0.01 m, while the considered dimensions for the straight beams' cross section are 0.01 m, 0.025 m and 0.05 m for both a and b.
The hexagonal grids are obtained by the assembly of 9 curved beams (see Fig. 2(a)).The overall dimensions of the hexagonal cell are kept constant.The only factor that is varied is the internal angle (β in Fig. 2(a)).The considered angles are β = 30 • , which corresponds to a regular hexagon, and β = −30 • .Cells with negative internal angles are considered to have re-entrant, or auxetic [13], geometry.The finite element configuration and geometry of both grids are summarized in Table 6.

Phase constant surfaces
The phase constants for the rectangular grids are first evaluated.A convergence study is performed to determine the required number of elements to obtain an accurate representation of the dynamic behavior of the cell in the considered frequency range.The results of this convergence study in shown in Fig. 7, which compares the phase constant estimated with a varying number of elements per side of the cell.In Fig. 7, the dashed line corresponds to a 1-element per side discretization, while the dotted and the thick solid line respectively show the results for 3 and 5 elements.Convergence is achieved by using 5 elements, as no changes in the phase constant configurations could be observed by refining the mesh beyond this value.Hence, in all subsequent analyses on the rectangular grid, a resolution of 5 elements is used as a good compromise between accuracy and computational efficiency.
The contour maps for the first phase surface for the considered rectangular grids are shown in Fig. 8.The isofrequency lines are labelled with the corresponding values of frequency expressed in rad/s.As discussed in Section 1, and in [5], the perpendicular to each isofrequency contour indicates the direction of wave propagation for the considered pair of propagation constants.Figure 8 shows that for frequencies below approximately 90 rad/s, the perpendicular to the iso-frequency contour lines spans an angular range between 0 and 90 0 as the propagation constants change.However, a transition occurs between 90 and 101 rad/s, above which the perpendicular to the iso-frequency lines spans a limited angular range.Accordingly, waves above the transition frequency tend to propagate only in limited directions and their influence is restricted to corresponding regions of the structure.This unique behavior of 2D periodic structures has already been demonstrated and discussed by Langley [3,4] and Ruzzene [5], and suggests possible interesting applications in systems or components where the confinement of waves and vibrations is essential to an effective operation.According to Langley and Ruzzene, the extension and in general the characteristics of the areas of confinement for elastic waves depend uniquely on the configuration of the unit cell in the 2D periodic structure.In here, the influence of various cross sectional areas for the beams along the circumferential direction is evaluated.Figure 8 shows that modifying the cross section of the curved beams significantly alters the phase constant surfaces.In the grid with a = 0.01 m, above the transition frequency ( ∼ = 90 rad/s), wave propagation mostly occurs along the circumferential direction (θ), while for a = 0.025 m and above the transition frequency ( ∼ = 210 rad/s), waves mostly propagate along the longitudinal direction (x).The iso-frequency lines above about 233 rad/s tend to be almost vertical and parallel to the θ axis, thus indicating strong directionality of wave motion along the longitudinal axis of the cylinder.This behavior is also observed at a lower frequency in the case where a = 0.050 m, where above transition ( ∼ = 100 rad/s), wave propagation occurs preferentially along the axial direction (x).

Harmonic response
The predictions obtained from the analysis of the phase constant surfaces are validated by evaluating the harmonic response of the complete 2D periodic assembly.The response of the first (a = 0.01 m) and third (a = 0.05 m) rectangular grid is considered for harmonic excitations of unit amplitude and frequency equal to 80, 105 and 125 rad/s.The load is applied at a location identified by the following set of cylindrical coordinates: θ = π, R = 1 m and z = 5 m.The results are shown in Figs 9-10, which show the magnitude of the out of plane motion of the grid through bi-dimensional maps, obtained by interpolating the nodal displacements resulting from the FE analysis.Figure 9(a) shows the harmonic response of the first grid at an excitation frequency of 80 rad/s.According to the phase constant surface shown in Fig. 9(a), at this frequency waves are free to propagate in both the circumferential and longitudinal

Performance of hexagonal grids 5.3.1. Phase constant surfaces
The directional pattern of the considered class of hexagonal grids can be modified by varying the topology of the unit cell, through changes in the value of the internal angle β (see Fig. 2).This is demonstrated by the phase constant surfaces shown in Fig. 11, which compare the characteristics of grids with β = +30 • and β = −30 • .In the grid with β = 30 • , above the transition frequency ( ∼ = 250 rad/s), wave propagation mostly occurs along the circumferential direction (θ), while for β = −30 • this transition occurs at a much lower frequency ( ∼ = 150 rad/s).In both cases the iso-frequency lines above transition tend to be almost horizontal and parallel to the ε x axis, thus indicating strong directionality of wave motion along the circumferential direction.
The convergence of the hexagonal phase constant surfaces are also investigated using the same procedure described in Section 5.2.It is found that the convergence occurs when 9 elements are used in each cell.The figure showing the results of the convergence has been omitted for the sake of brevity.

Harmonic response
The predictions obtained from the analysis of the phase constant surfaces are again verified by evaluating the harmonic response of the complete 2D periodic assembly.The response of both hexagonal grids (β = 30 • and β = −30 • ) are considered at frequencies below and above the transition.The phase constant surfaces analysis has indicated that for the first grid transition occurs at ∼ = 250 rad/s, while for the second grid it takes place at ∼ = 150 rad/s.The structure is again excited at a location along its center rim by a harmonic force of unit amplitude, and the corresponding displacement maps are shown in Figs 12 and 13. Figure 12(a) shows the harmonic response of the first grid to an excitation frequency of 75 rad/s, i.e. below transition, while the displacement map shown in Fig. 12(b) corresponds to an excitation frequency of 300 rad/s.At this frequency above transition, one can observe wave attenuation in the longitudinal direction, and propagation in the circumferential direction.Similar characteristics are found in the response of the second grid, which is here evaluated at 75 and 190 rad/s (see Fig. 13) The phenomenon of wave attenuation in the cylindrical grid is greatly affected by the geometry of the unit cell.By changing the internal angle of the hexagonal cell, the frequency at which directivity of wave propagation is reduced by ∼ = 100 rad/s despite keeping all other factors constant.

Conclusions
This work investigates the vibration and wave propagation characteristics of cylindrical grid-like structures.The combination of the theory of 2D periodic structures, rotationally periodic structures, and the FE method, gives a methodology by which the dynamic response of cylindrical grids can be effectively predicted.In particular, the  presented methodology allows the evaluation of location and extension of propagation and attenuation regions within the structure.The considered grids are composed of curved beam elements assembled over a cylindrical surface according to various topological configurations.The influence of the cell geometry and topology on the dynamics of wave propagation is investigated through the developed techniques.The presented results demonstrate the ability to modify the wave propagation characteristics of the grids by proper design of the unit cell.In addition, the results demonstrate that propagation in the circumferential or longitudinal directions can occur preferentially above a transition frequency value.
The methodology used for the study provide the guidelines for the design of stiffened cylindrical structures with enhanced attenuation capabilities.The application of the concepts presented in this work could therefore be extremely relevant for improving the dynamic performance of shell structures commonly used in aerospace or naval applications.In particular, design configurations for grid stiffened structure can be identified so that the transmission of vibrations towards specified locations and at certain frequencies is minimized.The study can be extended to include the routines which identify the geometry and topology of the unit cell needed to achieve desired transmissibility levels in specified directions and for given excitation frequencies.

 
where, Derivatives of v and w are obtained by substituting v and w for u.In Eq. ( 40), i = 1 for u, i = 2 for v, and i = 3 for w.Hence, if the above relation is applied to Eq. ( 40 where, h k The constitutive matrix, C, used for the computation of the mass and stiffness matrices of the beam element is as follows: where E is the Young's modulus and G is the shear modulus of the material.

Appendix B. Beduction matrices
The reduction matrices, A and B, are defined as follows: where I is the identity matrix.

Table 1
Properties of Cantilevered Curved Beam

Table 3
Material and geometry for the strip of 5 rectangular cells

Table 5
Geometry and material properties of cylinder