Mathematical Problems in Engineering
Volume 2016 (2016), Article ID 9651430, 15 pages
http://dx.doi.org/10.1155/2016/9651430

Research Article

Periodic Properties of 1D FE Discrete Models in High Frequency Dynamics

A. Żak, M. Krawczuk, and M. Palacz

Faculty of Electrical and Control Engineering, Gdańsk University of Technology, Ulice Narutowicza 11/12, 80-234 Gdansk, Poland

Received 16 September 2015; Revised 11 December 2015; Accepted 13 December 2015

Academic Editor: Gen Qi Xu

Copyright © 2016 A. Żak et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Finite element discrete models of various engineering 1D structures may be considered as structures of certain periodic characteristics. The source of this periodicity comes from the discontinuity of stress/strain field between the elements. This behaviour remains unnoticeable, when low frequency dynamics of these structures is investigated. At high frequency regimes, however, its influence may be strong enough to dominate calculated structural responses distorting or even falsifying them completely. In this paper, certain computational aspects of structural periodicity of 1D FE discrete models are discussed by the authors. In this discussion, the authors focus their attention on an exemplary problem of 1D rod modelled according to the elementary theory.

1. Introduction

Periodic structures due to their unique dynamic properties are widely applied in many real engineering structures, such as isolation devices, acoustic filters, and dampers. The dynamic behaviour of periodic structures, especially the existence of frequency band gaps and negative refraction indexes, has been under investigation by many researchers [1–5]. Many analytical and discrete methods have been developed and employed for the analysis of dynamic characteristics of periodic structures.

Based on the literature of the subject, the following list of the applied methods can be made: the lump mass method [6], the multiple scattering method [7], the transfer matrix method [8], the plane wave expansion (PWE) method [9], the finite difference method (FDM) [10], the finite element method (FEM) [11], the wavelet method [12], the Dirichlet-to-Neumann map method [13], the boundary element method (BEM) [14], as well as the spectral finite element method in the frequency domain (FD-SFEM) [15], and the spectral finite element method in the time domain (TD-SFEM) [16].

It should be realised that the methods mentioned above have numerous limitations. Analytical methods such as the lump mass method, multiple scattering method, and the FD-SFEM are not suitable for investigation of geometrically complex structures. These limitations are nonexistent in the case of discrete methods such as the FEM and the TD-SFEM. On the other hand, the results of numerical simulations obtained by the use of the FEM or the TD-SFEM are not as accurate as the results obtained by the use of analytical methods. Their accuracy can be improved by an adaptive increase in the number of FEs (-method) or in the degree of approximation polynomials (-method) [17]. It is worth noting that the accuracy of - and -methods has been studied extensively by many researchers in the context of global, local, pollution, or dispersion errors [18–22]. Interesting results on the comparison between adaptive - and -methods and nonuniform rational B-splines (NURBS) can be found in [23, 24]; however, periodic properties of the numerical models used have stayed beyond such interest.

It is known that an increase in the number of FEs (-method) is not as effective as an increase in the degree of approximation polynomials (-method). In the case of equidistant node, distribution within FEs a significant increase in the degree of approximation polynomials may lead to their uncontrollable oscillations near the edges of FEs, known as Runge’s phenomenon [25, 26]. For this reason in numerical simulations by use of the FEM, the degree of approximation polynomials is typically not higher than 3. Contrary to this, Runge’s phenomenon is not observed in the case of the TD-SFEM, which is based on nonequidistant node distributions. The degree of approximation polynomials used by the TD-SFME is theoretically unlimited, but in numerical simulations mainly Chebyshev or Legendre polynomials of the fifth degree are applied. A further increase in the degree of approximation polynomials improves the accuracy of numerical results, but visible frequency band gaps appear at the high end of the natural frequency spectra [26]. This effect is very similar to the behaviour of periodic structures, but the source of structural periodicity has a different nature.

The source of this periodicity comes from the discontinuity of stress/strain field between FEs [17, 27]. Thus, the silent assumption about the representation of FE discrete models of continuous media/structures is not valid, despite the fact that it is not always visible based on their structural dynamic responses at low frequency regimes. However, there are certain regimes (i.e., high frequency dynamics and wave propagation) where this silent assumption may lead to significant errors due to the fact that overlooked structural periodicity of FE discrete models strongly manifests itself, as presented in this paper.

In this paper certain computational aspects of structural periodicity of 1D FE discrete models are discussed by the authors. Discrete models characterised by equidistant node distribution and nonequidistant distributions based on the roots of Chebyshev and Legendre polynomials were investigated. As approximation functions, Chebyshev and Legendre as well as cubic B-spline polynomials were tested. In this discussion, the authors focus their attention on an exemplary problem of 1D rod modelled according to the elementary theory of rods.

2. Finite Element Models as Periodic Structures

Bloch’s theorem is a very powerful tool used by physicists to investigate the behaviour of electrons in crystalline structures. However, Bloch’s theorem is applicable in a much broader context. This is thanks to the dual nature of electrons that can be treated as matter waves, also known as de Broglie waves. In fact, Bloch’s theorem can be successfully applied not only to explain the behaviour of wave related phenomena in periodic media, such as photonic and phononic crystals, but also to explain the behaviour of other types of mechanical structures characterised by certain periodic properties.

In a 1D case, Bloch’s theorem states that solution to a wave propagation problem in a periodic medium, also known as a Bloch wave, can be split up into the product of two waves: periodic wave of the same periodicity as the medium and plane-wave [28]:

where denotes the wave number and is the imaginary unit.

A 1D bimaterial periodic medium built up from periodically spaced cells is presented in Figure 1. The lengths of particular material regions within a single cell can be denoted as and . These regions are characterised by different values of elastic moduli and as well as different values of material densities and . Thus, these two regions are also characterised by different values of wave propagation phase velocities and or, alternatively, by different values of wave numbers and , but they are associated with common frequency through the relation:

where phase velocities and are expressed as

Figure 1: A schematic representation of a 1D bimaterial periodic medium.

For cells spaced at regular intervals the periodicity of the medium requires that and therefore it can be written as

Despite the fact that Bloch’s theorem assumes infinite periodic media, it can also provide valuable information about the behaviour of finite structures built up from a number of periodically spaced cells, assuming that . This can be achieved under assumption of cyclic boundary conditions in a similar manner as in the case of Fourier series expansion of aperiodic functions. It can be written as

leading to the following relation:

where denotes the total number of cells within the structure under consideration. In practice at least 100 cells are required in order to reveal the periodic nature of structural characteristics without the assumption about boundary conditions.

The characteristic equation of the problem, binding together cell properties and their periodic distribution through wave numbers , , and , can be obtained in a relatively simple manner after necessary mathematical manipulation. For this, it should be assumed that both Bloch’s waves and and their derivatives remain continuous across the boundaries of the cells, as well as associated periodic waves and across the boundaries of the two materials [28]:

where periodic waves and can be expressed as follows:

This leads to the characteristic equation of the following form:

which defines certain bands of frequencies , within which real solutions exist or, alternatively, within which wave propagation phenomena can take place.

In the case of periodic media made out of a single material characterised by discontinuity of its properties between the cells, (9) can be further simplified under assumption that and at , which can be written as

where is a certain nondimensional constant that sets the intensity of the discontinuity [28].

A typical frequency spectrum of a 1D periodic medium, calculated based on (9), is presented in Figure 2. It reveals a number of frequency band gaps within the spectrum that are denoted as , where no solutions exist. This behaviour is a direct consequence of the discontinuity of the wave propagation phase velocities between the cells, which may result from different values of elastic moduli and/or material densities. It should be noted that frequency band gaps appear at multiples of total number of cells .

Figure 2: A frequency spectrum of a 1D periodic medium (km/s, , m, and ).

At this point, it should be understood that FE discrete models of various engineering structures may be considered as structures of certain periodic characteristics. The source of the periodicity comes from the discontinuity of stress/strain field between FEs themselves [17, 26]. This behaviour remains unnoticeable at low frequency regimes, when low frequency dynamics of these structures is investigated. At high frequency regimes, however, its influence may become strong enough to dominate calculated structural responses distorting or even falsifying them completely.

The intensity of the stress/strain field discontinuities can be straightforwardly estimated for equidistant node distributions within FEs. In this case, the 1D medium under investigation can be assumed as modelled by FEs of equal length based on the elementary rod theory. In order to conform with (10), the periodic type of boundary conditions is used. Stress and Strain within a particular FE are linked together by the following very well-known formula:

where denotes the elastic modulus, is a FE solution representing the longitudinal displacement, and is a coordinate in the local coordinate system of the element, while is the total number of FEs of a discrete model. The symbol denotes the 1st derivative in respect of .

It can be seen from (11) that the stress/strain discontinuities between two adjacent FEs and connected together at node are in the current case related to the discontinuity of displacement derivative due to finite degree of approximation polynomials used. It should be noted that this problem is not present in the FE formulation based on the so-called analytical shape functions [29].

For different degrees of approximation polynomials , the discontinuity of the displacement derivative can be evaluated based on its Taylor series expansion around node . In the case of the degree of approximation polynomial , it can be written as

where all higher order terms are zero due to the assumed degree of approximation polynomial , while may be considered here as the distance between element nodes. Noting that and based on (12), it can be found that the discontinuity of displacement derivative at node can be expressed as

where the superscripts and have been omitted, where and is the central difference representation of the 2nd order derivative of exact solution of the problem under investigation supported on 3 nodes of 2 adjacent FEs, as shown in Figure 3.

Figure 3: Analytical solution and its FE representations and for the degree of approximation polynomial .

In a very similar manner, the discontinuity of displacement derivative at node for higher degrees of approximation polynomials can be obtained as

where is the central difference representation of the 4th order derivative of the exact solution supported on 5 nodes of 2 adjacent FEs, as well as

and where and are the central difference representations of the 6th and 4th order derivatives of exact solution supported on 7 nodes of 2 adjacent FEs, respectively.

The discontinuity of displacement derivative calculated above can be conveniently expressed in a nondimensional manner, which allows one to use it together with (10) in place of constant . In such a case, it may be expressed as

It should be realised that for higher degrees of approximation polynomials or for nonuniform distribution of nodes within Fes, corresponding relations for the discontinuity of displacement derivatives are much more complex and require additional computational effort.

Frequency band gaps of the 1D periodic medium under investigation can be calculated analytically, based on (10) and (16), or numerically by the FEM for various degrees of approximation polynomials .

The results on analytical calculations presented in Figure 4 clearly indicate that the frequency band gaps in the natural frequency spectra of FEM numerical models have their source in the stress/strain discontinuity between adjacent FEs. The number of frequency band gaps that can be observed is directly related to the degree of approximation polynomials and equal to . So there are no frequency band gaps () in the calculated frequency spectra for the FEM numerical model based on the first degree of approximation polynomials , while in the case of the third degree of approximation polynomials there are two frequency band gaps ().

Figure 4: Frequency band gaps of a 1D periodic medium, calculated analytically based on (10) and (16), for various degrees of approximation polynomials (km/s, m).

The widths of frequency band gaps are also directly related to the degree of approximation polynomials . They appear for wavelengths equal to the FE length or its multiples/fractions, that is, for at or for at . However, it should also be noted that for multiple frequency band gaps associated with degrees of approximation polynomials those band gaps that are related to long wavelengths are small and can be practically neglected—as seen from Figure 4 in the case of in comparison with .

The results presented in Figure 5 correspond to numerical calculations by the FEM for the same three degrees of approximation polynomials . They correspond very well to the results presented in Figure 4, but the frequency band gaps observed are broader than those predicted analytically based on exact solution . This is due to overstiffening of FEM solutions above the Nyquist frequency [25, 30] because above this frequency FEM solutions try to represented a continuous analytical solution that is based on infinite degrees of freedom by a discrete numerical solution based on the finite number of degrees of freedom.

Figure 5: Frequency band gaps of a 1D periodic medium, calculated numerically by the FEM, for various degrees of approximation polynomials —a numerical model of 300 DOFs according to the elementary theory of rods [36], based on equidistant node distribution (km/s, m).

It turns out that the behaviour observed in Figures 4 and 5 stays closely correlated not only with the degree of approximation polynomials , but also with the element formulation (stress or displacement). These aspects of FE modelling are discussed in more detail in the following sections of this paper. As before, frequency band gaps appear at multiples of total number of cells/elements within the discrete numerical model.

3. Numerical Simulations

All results of numerical calculations presented in the following parts of this paper were selected by the authors in order to reveal the periodicity nature of FE models in the most legible manner. However, it should be noted that in a general case FE models employed nowadays to solve numerous scientific problems combine together various types of FEs (1D, 2D, or 3D) and material properties (isotropic, orthotropic, or anisotropic), as well as discretisation techniques (regular, spectral, or adaptive). This results in a very complex and coupled numerical behaviour, in which particular influences related to either modelling or discretisation are difficult to separate from each other. In fact a great variety of FE techniques used [17, 26, 29, 31, 32] (time domain, frequency domain, classical, spectral, or semianalytical) additionally adds up to the complexity of the problem under consideration.

In order to avoid such complex and coupled effects, the authors decided to focus attention on an exemplary 1D problem. In this problem, an isotropic rod is considered that is modelled according to the elementary theory of rods [29, 33]. Since this rod theory results in a nondispersive behaviour [29, 33], thanks to this any dispersion effects observed in results of numerical simulations must have numerical origins.

The geometry of the rod under investigation is shown in Figure 6. It is assumed that the rod is made out of aluminium alloy of elastic modulus GPa, material density kg/m³, and remains free—its both ends remain unsupported. The length of the rod is m, while its height and width are the same mm. For the assumed rod theory, phase and group velocity of propagating elastic waves are the same. For the given material properties, they are equal to km/s. In such a case, longitudinal waves propagating within the rod require 1 ms to travel over length of the rod from point P₁ to point P₂.

Figure 6: Geometry of an isotropic rod of free ends.

At first, the analysis of rod natural frequency spectra for various degrees of approximation polynomials was performed. The analysis aimed to reveal the periodic behaviour of FE models related to their discretisation into FEs. In the analysis, numerical models of 736 DOFs were used for Chebyshev node distribution within the elements. It can be clearly seen from Figure 7 that an increase in the degree of approximation polynomials increases the number of frequency band gaps, denoted as . Their number is directly related to the degree of approximation polynomials and is equal to . As a consequence, for linear approximation polynomials no frequency band gaps are observed. Moreover, the frequency band gaps in the part of the frequency spectra corresponding to lower natural frequencies are relatively small, therefore in most cases their presence may remain unnoticed.

Figure 7: Frequency spectra of an isotropic rod of free ends for various degrees of approximation polynomials —a numerical model of 736 DOFs according to the elementary theory of rods [36], based on Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

Based on the results presented in Figure 7, it may seem reasonable to reduce the degree of approximation polynomials down to to get rid of any unwanted periodic behaviour of FE models. It is well known [22] that in such a case the approximation error is strongly increased in comparison with higher degrees of approximation polynomials . Additionally, it should be noticed that the degree of approximation polynomials corresponds to equidistant node distribution characteristic for the classical FEM.

This behaviour is well illustrated by Figure 8, in which relative errors of rod natural frequency spectra are presented for the same degrees of approximation polynomials , which is the second part of the analysis. Relative error was calculated based on well-known analytical formula [34] for the type of boundary condition considered:

where denote values of natural frequencies calculated numerically, while is their corresponding analytical values.

Figure 8: Relative errors of frequency spectra of an isotropic rod of free ends for various degrees of approximation polynomials —a numerical model of 736 DOFs according to the elementary theory of rods [36], based on Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

It is well seen from Figure 8 that the use of higher degrees of approximation polynomials results in a significant increase in relative error in the part of the frequency spectra corresponding to higher natural frequencies, where frequency gaps become considerable. Despite this fact in the part of the frequency spectra corresponding to lower natural frequencies, FE models based on higher degrees of approximation polynomials are characterised by much lower values of relative error .

This error depends not only on the degrees of approximation polynomial , but also on the node distributions (nonequidistant or equidistant) within FEs associated with the type of approximation polynomials (nonequidistant Chebyshev and Legendre and equidistant B-spline), as presented in Table 1. Errors and were defined in a very simple manner as

Table 1: Average errors of frequency spectra of an isotropic rod of free ends for various types and degrees of approximation polynomials —a numerical model of 736 DOFs according to the elementary theory of rods [36] used for calculations by the TD-FEM ( km/s,  m).

Polynomial type
Polynomial degree
[%]
[%]

Chebyshev
1
3.47
9.92

3
0.44
11.73

5
0.13
15.36

7
0.05
18.75

Legendre
1
3.39
12.76

3
0.25
3.81

5
0.06
8.95

7
0.02
13.01

B-spline
3
0.01
0.85

The results collected in Table 1 indicate that the type of approximation polynomials plays a significant role as a factor determining the accuracy of FE models. In the case of nonequidistant distributions of nodes within elements, based on either Chebyshev or Legendre polynomials [26], it can be noted that for the same degrees of approximation polynomials errors and are smaller for approximation polynomials based on Legendre node distribution. As before, it should be emphasised that the degree of approximation polynomials corresponds to equidistant node distribution characteristic for the classical FEM.

In the part of frequency spectra considered for computational practice and corresponding to lower natural frequencies, error decreases rapidly with an increase in the degree of approximation polynomial . In the case of approximation polynomials based on Chebyshev node distribution, error decreases 73 times, while in the case of Legendre polynomials it decreases 192 times. Lower errors and observed for Legendre node distribution are a direct consequence of the numerical integration scheme used for computation of elemental characteristic inertia matrices by Lobatto quadrature that leads to their diagonal form. Lobatto quadrature underestimates values of elemental characteristic inertia matrices [35] as exact up to the degree of approximation polynomials equal to .

Gauss quadrature employed in the case of Chebyshev node distribution is exact up to the degrees of approximation polynomials equal to but leads to full forms of elemental characteristic inertia matrices. For the same quadrature B-spline cubic approximation, polynomials based on equidistant node distribution evidently offer the best results. This is because of the fact that B-spline approximation polynomials enforce not only the continuity of approximated functions, but also the continuity of their derivatives up to degree, which results in no frequency gaps observed since stress/strain fields between FEs remain continuous. The natural frequency spectrum of the rod under consideration for B-spline cubic approximation polynomials and equidistant node distribution is presented in Figure 9, while corresponding relative error is presented in Figure 10.

Figure 9: Frequency spectrum of an isotropic rod of free ends for B-spline cubic approximation polynomials—a numerical model of 736 DOFs according to the elementary theory of rods [36], based on equidistant node distribution used for calculations by the FEM (km/s, m).

Figure 10: Relative error of the frequency spectrum of an isotropic rod of free ends for B-spline cubic approximation polynomials—a numerical model of 736 DOFs according to the elementary theory of rods [36] based on equidistant node distribution used for calculations by the FEM (km/s, m).

4. Known Issues and Their Consequences

The third part of the current analysis is concerned with certain numerical issues that are typical examples of the problems related with the application of the FEM or the TD-SFEM. In this analysis, force vibration responses and time responses of the rod under investigation were also examined.

Two types of rod discretisation were selected for the analysis. In the first case, a numerical model of the rod was based on Chebyshev node distribution and approximation polynomials of the fifth degree. In the second case, the corresponding numerical model of the rod was based on equidistant node distribution and B-spline cubic approximation polynomials. As before a numerical model of 736 DOFs according to the elementary theory of rods was employed.

In order to conform with the results presented in Figures 7 and 9, force vibration responses of the rod were calculated assuming an excitation force variable within a frequency band from 0 Hz to 350 kHz and acting at point P₁  (), where also the amplitude of rod vibration was assessed. The amplitude of the excitation was 1 N. A small value of material damping was assumed, expressed by the loss factor , to bound vibration responses of the rod nearby its resonant frequencies.

Within the assumed frequency band frequency band gap FG₃, around 230 kHz should most prominently affect rod dynamic behaviour. For this reason, the results of numerical calculations by the FEM presented in Figure 11 were narrowed to a frequency band from 200 kHz to 250 kHz.

Figure 11: Forced vibration response of an isotropic rod of free ends within a frequency band from 200 kHz to 250 kHz for Chebyshev approximation polynomials—a numerical model of 736 DOFs according to the elementary theory of rods [36], based on Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

It is evident that the presence of frequency band gap FG₃ results in inability of the rod to produce any dynamic responses within a 13 kHz frequency band gap affected by the gap and starting from 223.5 kHz up to 236.5 kHz. Obviously such a feature of the rod FE model must strongly influence calculated time responses that are in fact built as sums of harmonic responses of various amplitudes. At this point, it is interesting to note that in the case of the rod numerical model based on equidistant node distribution and B-spline cubic approximation polynomials dynamic responses of the rod reveal no frequency band gaps due to the continuity of the stress/strain field between FEs, as presented in Figure 12.

Figure 12: Forced vibration response of an isotropic rod of free ends within a frequency band from 200 kHz to 250 kHz for B-spline cubic approximation polynomials—a numerical model of 736 DOFs according to the elementary theory of rods [36], based on equidistant node distribution used for calculations by the FEM (km/s, m).

During the analysis of time responses of the rod, various carrier frequencies of the excitation force acting at point P₁  were examined with their spectra gradually closing to frequency gap FG₃ [26]. As the excitation, a 12-cycle sine signal modulated by the Hann window was selected, while the excitation amplitude was 1 N. As before time responses of the rod were assessed at the excitation point. The total calculation time was 2.5 ms and was divided into as many as 25,000 equal time steps to minimise the time discretisation error. Numerical solution of the equation of motion was obtained by the application of the central difference method.

For approximation polynomials of the fifth degree based on Chebyshev node distribution, calculated time responses of the rod are presented in Figure 13. In this figure, times and represent analytically calculated time windows, within which the time responses should be visible under assumption of no numerical dispersion of the FE model employed. The times and are dependent on the excitation carrier frequency and their values are shown in Table 2.

Table 2: Spans of propagating signal time windows for various excitation carrier frequencies calculated analytically according to the elementary theory of rods [36] ( km/s,  m).

Carrier frequency [kHz]
[ms]
[ms]

75
2.0
2.160

100
2.0
2.120

125
2.0
2.096

150
2.0
2.080

Figure 13: Time responses of an isotropic rod of free ends for various excitation carrier frequencies for approximation polynomials of the fifth degree—a numerical model of 736 DOFs according to the elementary theory of rods [36], based on Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

The results obtained clearly indicate that for excitation carrier frequencies distance from frequency band gap FG₃, up to 100 kHz, the observed wave propagation patterns remain undistorted. For higher excitation carrier frequencies, gradually approaching the frequency band affected by the presence of frequency gap FG₃, these patterns become gradually deformed. In the case of the exaction carrier frequency equal to 150 kHz, the misrepresentation of the wave propagation patter obtained is very clearly visible.

The reason for such behaviour comes from the fact that time response frequency components close to frequency band gap FG₃ are strongly affected by an increasing numerical modelling error, as shown in Figure 8. For example, the frequency content of 150 kHz excitation is presented in Figure 14.

Figure 14: Frequency spectrum of 150 kHz excitation in the form of 12-cycle sine signal modulated by the Hann window.

As long as the excitation frequency components fall into the frequency band influenced by frequency gaps, calculated time responses of the rod must carry corresponding information that manifests itself as numerical dispersion. Certainly the intensity of this numerical dispersion is very small in the case of excitation frequencies from the lower part of frequency spectra and therefore can be neglected. In the case of excitation frequencies from the upper part of frequency spectra, numerical dispersion can be significant and in extreme cases, in the close vicinity of frequency band gaps, may prevent obtaining any numerical solution.

Contrary to this in the case of the rod numerical model based on equidistant node distribution and B-spline cubic approximation, the numerical dispersion problem described above is not present due to no frequency gaps observed in the rod natural frequency spectrum. This is illustrated by Figure 15. As a consequence of that, no numerical dispersion resulting from rod modelling is observed and calculated wave propagation patterns remain undistorted.

Figure 15: Time responses of an isotropic rod of free ends for various excitation carrier frequencies for B-spline cubic approximation polynomials—a numerical model of 736 DOFs according to the elementary theory of rods [36] based on equidistant node distribution used for calculations by the FEM (km/s, m).

Certainly, it can be expected that the manifestation of the periodicity in the case of complex FE discrete models combining together various types of FEs (1D, 2D, or 3D), material properties (isotropic, orthotropic, or anisotropic), as well as discretisation techniques (regular, spectral, or adaptive), must result in a more complex behaviour. However, it can be also expected that certain periodic features of this behaviour should remain the same. This can be well illustrated based on the result presented in Figures 16 and 17 that present frequency spectra of the same rod obtained for the two-mode Mindlin-Herrmann theory of rods. The symbols S₀ and S₁ denote the first and second symmetrical (longitudinal) wave propagation modes, while and are in fact material constants expressing the velocities of longitudinal and shear waves in 3D elastic media, respectively [36].

Figure 16: A frequency spectrum of an isotropic rod of free ends for approximation polynomials of the fifth degree—a numerical model of 1502 DOFs according to the two-mode Mindlin-Herrmann theory of rods [36], based on Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

Figure 17: A frequency spectrum of an isotropic rod of free ends for B-spline cubic approximation polynomials—a numerical model of 1502 DOFs according to the two-mode Mindlin-Herrmann theory of rods [36], based on equidistant node distribution used for calculations by the FEM (km/s, m).

Two different discretisation methods of the rod were taken into account. In the first case, the fifth degree of approximation polynomials for Chebyshev node distribution was employed, while in the second case they were B-spline cubic approximation polynomials for equidistant node distribution. Both models included 1502 DOFs.

Figures 16 and 17 show that the two-mode Mindlin-Herrmann theory requires that the frequency spectrum of the rod must be split between these two modes. Thus, above a certain frequency, known as the first cut-off frequency , equal to 419.3 kHz in the current case, both the modes are simultaneous—please refer to Figure 18 for more details. Thus, those rod natural frequencies, the same as natural vibration modes, are dual above the cut-off frequency .

Figure 18: Dispersion curves for the group velocity for the two-mode Mindlin-Herrmann theory of rod (km/s, km/s).

It can be seen that in the case of the rod FE discrete model according to the two-mode Mindlin-Herrmann theory, based on the fifth degree of approximation polynomials and Chebyshev node distribution, due to the discontinuity of the stress/strain field between FEs previously discusses, the number of frequency band gaps is also doubled in comparison with the elementary rod theory.

This has profound consequences on the accuracy of numerical simulations, since the spans of particular frequency band gaps increase according the wave propagation mode and not according to their absolute number—the span of frequency band gap FG₄ just below the first cut-off frequency is much bigger than frequency band gap FG₅. As before no frequency band gaps are observed for B-spline cubic approximation polynomials and equidistant node distribution. In this case the separation of the frequency spectrum between the two modes below and above the first cut-off frequency is very clear.

Finally, two important issues related to FE discrete numerical models should be addressed. The results of numerical calculations presented in Figure 19 clearly indicate that the observed behaviour related to the presence of frequency band gaps in the frequency spectrum of the 1D rod under investigation should be attributed to the periodic properties of its FE discrete model. Moreover, this model feature should not be understood as a numerical error but a numerical feature. It can be seen from Figure 19 that regardless the size of rod numerical model, in terms of its elements/DOFs, the same general pattern is observed. The presence of frequency band gaps in the frequency spectrum is directly linked to the degree of approximation polynomials employed. Moreover, the frequency spectrum is effectively divided into segments with frequency band gaps between them.

Figure 19: Relative errors of frequency spectra of an isotropic rod of free ends for various numbers of FEs used for rod discretisation—a numerical model according to the elementary theory of rods [36], based on Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

In the low frequency regime, no periodic features in the rod frequency spectrum are seen. For this reason, this part of the spectrum may be considered as representative for the dynamic behaviour of continuous media/structures. Contrary to that the upper part of the spectrum must be considered as affected by the periodic properties of the numerical discrete model itself. This becomes a very important matter when periodic structures are studied, modelled by the use of the FEM or the TD-SFEM.

The second issue comes from the obvious fact that a majority of commercial FE packages offers automatic mesh generation capabilities in the case of 1D, 2D, and 3D problems. Such automatic mesh generators produce irregular meshes that are characterised by some distribution of element dimensions, which are FE lengths in the case of 1D problems under investigation.

However, it should be understood that the periodic properties of FE discrete models have their origins in the stress/strain discontinuities between adjacent FEs. As long as the dimensions of FEs concentrate around a common size, which can be understood as the average element size, one can expect that certain characteristics that can be attributed to the periodicity of FE discrete models should appear. This is well illustrated by Figure 20.

Figure 20: Relative errors of frequency spectra of an isotropic rod of free ends for various dispersions of element length —a numerical model of 736 DOFs according to the elementary theory of rods [36], based on the fifth degree of approximation polynomials and Chebyshev node distribution used for calculations by the TD-FEM (km/s, m).

It was assumed that the lengths of FEs used for modelling the rod under consideration could vary within a certain limit, denoted as , but still these lengths are concentrated around a common average value . It can be noted from Figure 20 that an increase in the dispersion of rod element lengths slightly affects the observed error of rod natural frequencies. The strongest influence is observed in the case of the dispersion of rod element lengths equal to . In this case, the frequency band gaps FG₁, FG₂, and FG₃ associated with the lower part of the rod natural frequency spectra are smoothed, but the high frequency behaviour remains the same; that is, the frequency band gap FG₄ stays visible; however, its frequency span is slightly smaller.

5. Conclusions

In this paper certain computational aspects of the periodic nature of FE discrete models have been investigated and discussed by the authors. The results of this analysis allow them to formulate certain general conclusions:

(i)The majority of FE discrete models of engineering structures used in the analysis of high frequency dynamics or wave propagation problems represent structures of certain periodic characteristics.(ii)This analysis requires uniform or nearly uniform meshes of FEs in order to avoid problems with artificial anisotropy of FE discrete models, typically assuming at least 5 element nodes per wavelength associated with the slowest propagating wave associated.(iii)The source of the periodicity comes from the fact that the classical formulation of the FEM or the TD-SFEM requires only the continuity of the displacement field between FEs, while the stress/strain field between the elements may remain discontinuous but finite.(iv)The discontinuity of the stress/strain filed between FEs results in the discontinuity of frequency spectra and the presence of the so-called frequency band gaps, which represent frequency domains, where no real solutions exist.(v)The number of frequency band gaps observed in the frequency spectra of discrete FE models is directly related to the number of FEs used as well as the degree of approximation polynomials employed.(vi)At low frequency regimes, when low frequency dynamic or static structural behaviour is investigated, the presence of frequency band gaps may remain unnoticed.(vii)At high frequency regimes, the presence of frequency band gaps may dominate calculated high frequency structural dynamic responses distorting or even falsifying them completely, especially in the case when multimode theories are employed.(viii)In the case of the classic formulation of the FEM or the TD-SFEM, lower degrees of approximation polynomials, regardless of available node distributions such as equidistant or nonequidistant Legendre or Chebyshev, result in a smaller number of frequency band gaps, while FE discrete models based on such approximation polynomials are characterised by much lower accuracy.(ix)Contrary to that, higher degrees of approximation polynomials result in a greater number of frequency band gaps, while corresponding FE discrete models based on such approximation polynomials are characterised by much higher accuracy.(x)In the case of higher degrees of approximation polynomials, the presence of frequency band gaps reduces the usable part of frequency spectra to their lower parts.(xi)For FE discrete models based on multimode theories, their periodic nature has the same origin and character; however, the observed behaviour is more complex.(xii)The influence of the periodicity on FE discrete models may be minimised or fully eliminated in the case of the stress (compliance) formulation of the FEM or in the case of the superelement approach of the FEM or the TD-SFEM.(xiii)The stress (compliance) formulation of FE discrete models, enforcing the continuity of the stress field between FEs, leads to the formulation of the characteristic elemental matrices equivalent to the case, when linear approximation polynomials are used, thus characterised by relatively low accuracy for higher frequency regimes.(xiv)The superelement approach, enforcing the continuity of stress/strain field within the whole superelement, becomes numerically inefficient for large discrete models as leading to the full stiffness matrix of the superelement.(xv)The formulation of the FEM based on B-spline cubic approximation polynomials, regardless of node distribution, fully eliminates the influence of the periodicity of FE discrete models due to the continuity of the stress/strain field between the elements and is free of the drawbacks associated with the stress or superelement formulations of the FEM or TD-SFEM.(xvi)Appropriate modelling of periodic structures by discrete FE models requires special attention in order to control and to separate the influence of periodic properties of the models from the periodicity of the structures themselves.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors of this work would like to gratefully acknowledge the support for their research provided by the National Science Centre through via Project UMO-2012/07/B/ST8/03741 Wave propagation in periodic structures. All results presented in this paper have been obtained by the use of the software available at the Academic Computer Centre in Gdańsk in the frame of a computational project.

References

M. S. Kushwaha, P. Halevi, L. Dobrzynski, and B. Djafari-Rouhani, “Acoustic band structure of periodic elastic composites,” Physical Review Letters, vol. 71, no. 13, pp. 2022–2025, 1993.

X. D. Zhang and Z. Y. Liu, “Negative refraction of acoustic waves in two-dimensional phononic crystals,” Applied Physics Letters, vol. 85, no. 2, pp. 341–343, 2004.

R. X. Feng and K. X. Liu, “Tuning the band-gap of phononic crystals with an initial stress,” Physica B: Condensed Matter, vol. 407, no. 12, pp. 2032–2036, 2012.

R. K. Narisetti, M. J. Leamy, and M. Ruzzene, “A perturbation approach for predicting wave propagation in one-dimensional nonlinear periodic structures,” Journal of Vibration and Acoustics, vol. 132, no. 3, Article ID 031001, 2010.

M. M. Sigalas and E. N. Economou, “Elastic and acoustic wave band structure,” Journal of Sound and Vibration, vol. 158, no. 2, pp. 377–382, 1992.

G. Wang, J. Wen, and X. Wen, “Quasi-one-dimensional phononic crystals studied using the improved lumped-mass method: application to locally resonant beams with flexural wave band gap,” Physical Review B, vol. 71, no. 10, Article ID 104302, 2005.

I. E. Psarobas, N. Stefanou, and A. Modinos, “Scattering of elastic waves by periodic arrays of spherical bodies,” Physical Review B, vol. 62, article 278, 2000.

C. Cetinkaya and A. F. Vakakis, “Transient axisymmetric stress wave propagation in weakly coupled layered structures,” Journal of Sound and Vibration, vol. 194, no. 3, pp. 389–416, 1996.

M. M. Sigalas and E. N. Economou, “Elastic waves in plates with periodically placed inclusions,” Journal of Applied Physics, vol. 75, no. 6, pp. 2845–2850, 1994.

M. M. Sigalas and N. García, “Theoretical study of three dimensional elastic band gaps with the finite-difference time-domain method,” Journal of Applied Physics, vol. 87, no. 6, pp. 3122–3125, 2000.

Y. Liu and L.-T. Gao, “Explicit dynamic finite element method for band-structure calculations of 2D phononic crystals,” Solid State Communications, vol. 144, no. 3-4, pp. 89–93, 2007.

Z.-Z. Yan and Y.-S. Wang, “Calculation of band structures for surface waves in two-dimensional phononic crystals with a wavelet-based method,” Physical Review B, vol. 78, no. 9, Article ID 094306, 2008.

F.-L. Li and Y.-S. Wang, “Application of dirichlet-to-neumann map to calculation of band gaps for scalar waves in two-dimensional phononic crystals,” Acta Acustica united with Acustica, vol. 97, no. 2, pp. 284–290, 2011.

F.-L. Li, Y.-S. Wang, and C. Zhang, “Boundary element method for bandgap computation of photonic crystals,” Optics Communications, vol. 285, no. 5, pp. 527–532, 2012.

Z.-J. Wu, F.-M. Li, and Y.-Z. Wang, “Study on vibration characteristics in periodic plate structures using the spectral element method,” Acta Mechanica, vol. 224, no. 5, pp. 1089–1101, 2013.

A. Żak, M. Krawczuk, and L. Doliński, “Dynamics of 1-D periodic structures. Numerical and experimental investigation,” in Proceedings of the 21st International Congress on Sound and Vibration, pp. 1–8, Beijing, China, July 2014.

O. C. Zienkiewicz, The Finite Element Method, Butterworth-Heinemann, Oxford, UK, 2000.

F. Ihlenburg and I. Babuška, “Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation,” International Journal for Numerical Methods in Engineering, vol. 38, no. 22, pp. 3745–3774, 1995.

J. T. Oden and Y. Feng, “Local and pollution error estimation for finite element approximations of elliptic boundary value problems,” Journal of Computational and Applied Mathematics, vol. 74, no. 1-2, pp. 245–293, 1996.

I. Babuška, T. Strouboulis, S. K. Gangaraj, and C. S. Upadhyay, “Pollution error in the h-version of the finite element method and the local quality of the recovered derivatives,” Computer Methods in Applied Mechanics and Engineering, vol. 140, no. 1-2, pp. 1–37, 1997.

A. Deraemaeker, I. Babuška, and P. Bouillard, “Dispersion and pollution of the FEM solution for the helmholtz equation in one, two and three dimensions,” International Journal for Numerical Methods in Engineering, vol. 46, no. 4, pp. 471–499, 1999.

M. Jasiński and G. Zboiński, “On some hp-adaptive finite element method for natural vibrations,” Computers & Mathematics with Applications, vol. 66, no. 11, pp. 2376–2399, 2013.

J. A. Cottrell, A. Reali, Y. Bazilevs, and T. J. R. Hughes, “Isogeometric analysis of structural vibrations,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 41–43, pp. 5257–5296, 2006.

T. J. R. Hughes, A. Reali, and G. Sangalli, “Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of p-method finite elements with k-method NURBS,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 49-50, pp. 4104–4124, 2008.

A. Zak and M. Krawczuk, “Certain numerical issues of wave propagation modelling in rods by the Spectral Finite Element Method,” Finite Elements in Analysis and Design, vol. 47, no. 9, pp. 1036–1046, 2011.

W. Ostachowicz, P. Kudela, M. Krawczuk, and A. Żak, Guided Waves in Structures for SHM. The Time-domain Spectral Element Method, John Wiley & Sons, Singapore, 2012.

J. N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill Education, New York, NY, USA, 2005.

D. J. Griffths, Introduction to Quantum Mechanics, Pearson Education, Harlow, UK, 2004.

J. F. Doyle, Wave Propagation in Structures, Springer, New York, NY, USA, 1997.

Y. M. Xie and G. P. Steven, “Explicit formulas for correcting finite-element predictions of natural frequencies,” Communications in Numerical Methods in Engineering, vol. 9, no. 8, pp. 671–680, 1993.

S. S. Rao, The Finite Element Method in Engineering, Pergamon Press, Oxford, UK, 1981.

Y. K. Cheung, Finite Strip Method in Structural Analysis, Pergamon Press, Oxford, UK, 1976.

J. W. S. Rayleigh, The Theory of Sound, Dover, New York, NY, USA, 1945.

R. D. Blevins, Formulas for Natural Frequency and Mode Shapes, Krieger Publishing, Amsterdam, The Netherlands, 1976.

A. Ralston, A First Course in Numerical Analysis, McGraw-Hill, New York, NY, USA, 1965.

A. Żak and M. Krawczuk, “Assessment of rod behaviour theories used in spectral finite element modelling,” Journal of Sound and Vibration, vol. 329, no. 11, pp. 2099–2113, 2010.