#### 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

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 asleading 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 aswhere 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 .

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 aswhere 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 aswhere 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.

In a very similar manner, the discontinuity of displacement derivative at node for higher degrees of approximation polynomials can be obtained aswhere is the central difference representation of the 4th order derivative of the exact solution supported on 5 nodes of 2 adjacent FEs, as well asand 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 ().

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.

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^{3}, 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_{1} to point P_{2}.

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.

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.

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

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.

#### 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_{1} (), 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_{3}, 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.

It is evident that the presence of frequency band gap FG_{3} 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.

During the analysis of time responses of the rod, various carrier frequencies of the excitation force acting at point P_{1} were examined with their spectra gradually closing to frequency gap FG_{3} [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.

The results obtained clearly indicate that for excitation carrier frequencies distance from frequency band gap FG_{3}, 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_{3}, 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_{3} 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.

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.

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_{0} and S_{1} 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].

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 .

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_{4} just below the first cut-off frequency is much bigger than frequency band gap FG_{5}. 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.

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.

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_{1}, FG_{2}, and FG_{3} 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_{4} 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.