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

Authors: A. Żak; M. Krawczuk; M. Palacz
Journal: Mathematical Problems in Engineering
 (2016)

---

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


---

## Body

## 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 (h-method) or in the degree of approximation polynomials (p-method) [17]. It is worth noting that the accuracy of h- and p-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 h- and p-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 (h-method) is not as effective as an increase in the degree of approximation polynomials (p-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 solutionf(x) 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 g(x) of the same periodicity as the medium and plane-wave eikx [28]:(1)fx=eikxgx,where k denotes the wave number and i is the imaginary unit.A 1D bimaterial periodic medium built up from periodically spaced cells is presented in Figure1. The lengths of particular material regions within a single cell can be denoted as (a-b) and b. These regions are characterised by different values of elastic moduli E1 and E2 as well as different values of material densities ρ1 and ρ2. Thus, these two regions are also characterised by different values of wave propagation phase velocities c1 and c2 or, alternatively, by different values of wave numbers k1 and k2, but they are associated with common frequency ω through the relation:(2)ω=k1c1=k2c2,where phase velocities c1 and c2 are expressed as(3)c1=E1ρ1,c2=E2ρ2.Figure 1
A schematic representation of a 1D bimaterial periodic medium.For cells spaced at regular intervalsa the periodicity of the medium requires that g(x+a)=g(x) and therefore it can be written as(4)fx+a=eikx+agx+a=eikaeikxgx=eikafx.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 ofN periodically spaced cells, assuming that N≫1. This can be achieved under assumption of cyclic boundary conditions f(x+aN)=f(x) in a similar manner as in the case of Fourier series expansion of aperiodic functions. It can be written as(5)fx+aN=fx=eikaNfx⟶eikaN=1leading to the following relation:(6)eikaN=1⟶k=2πnaN,n=1,2,…,N,where N 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 numbersk, k1, and k2, can be obtained in a relatively simple manner after necessary mathematical manipulation. For this, it should be assumed that both Bloch’s waves f1(x) and f2(x) and their derivatives remain continuous across the boundaries of the cells, as well as associated periodic waves g1(x) and g2(x) across the boundaries of the two materials [28]:(7)f10+=f20-,g1a-b=g2b,f1′0+=f2′0-,g1′a-b=g2′b,where periodic waves g1(x) and g2(x) can be expressed as follows:(8)g1x=A1eik1-kx+A2e-ik1+kx,g2x=B1eik2-kx+B2e-ik2+kx.This leads to the characteristic equation of the following form:(9)coska=cosk1a-bcosk2b-k12+k222k1k2sink1a-bsink2bwhich defines certainbands 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 b→0 and k2→∞ at k2b=const, which can be written as(10)coska=cosk1a-Asin⁡k1ak1a,where A 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 FG1,…,FG4, 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 FG1,…,FG4 appear at multiples of total number of cells N.Figure 2
A frequency spectrum of a 1D periodic medium (c1=5km/s, c2=0.5c1, a=0.025m, and b=0.05a).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 lengtha based on the elementary rod theory. In order to conform with (10), the periodic type of boundary conditions is used. Stress σj(x) and Strain ϵj(x) within a particular FE j are linked together by the following very well-known formula:(11)σjx=Eϵjx=Eu′jx,j=1,2,…,J,where E denotes the elastic modulus, uj(x) is a FE solution representing the longitudinal displacement, and x is a coordinate in the local coordinate system of the element, while J is the total number of FEs of a discrete model. The symbol ′ denotes the 1st derivative in respect of x.It can be seen from (11) that the stress/strain discontinuities between two adjacent FEs j and j+1 connected together at node s are in the current case related to the discontinuity of displacement derivative u′j∣j+1(xs) due to finite degree p 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 polynomialsp, the discontinuity of the displacement derivative u′j∣j+1(xs) can be evaluated based on its Taylor series expansion around node s. In the case of the degree of approximation polynomial p=1, it can be written as(12)ujxs-h=ujxs-u′jxs-huj+1xs+h=uj+1xs+u′j+1xs+h,j=1,2,…,J,where all higher order terms are zero due to the assumed degree of approximation polynomial p=1, while h may be considered here as the distance between element nodes. Noting that uj(xs)=uj+1(xs) and based on (12), it can be found that the discontinuity of displacement derivative u′j∣j+1(xs) at node s can be expressed as(13)Δsp=1=u′xs--u′xs+=-1hus-1-2us+us+1≈-fs∣32h,where the superscripts j and j+1 have been omitted, where us=u(xs) and fs∣3(2) is the central difference representation of the 2nd order derivative of exact solution f(x) of the problem under investigation supported on 3 nodes of 2 adjacent FEs, as shown in Figure 3.Figure 3
Analytical solutionf(x) and its FE representations uj(x) and uj+1(x) for the degree of approximation polynomial p=1.In a very similar manner, the discontinuity of displacement derivativeu′j∣j+1(xs) at node s for higher degrees of approximation polynomials p can be obtained as(14)p=2⟶Δsp=2≈12fs∣54h3,where fs∣5(4) is the central difference representation of the 4th order derivative of the exact solution f(x) supported on 5 nodes of 2 adjacent FEs, as well as(15)p=3⟶Δsp=3≈-512fs∣76h5+12fs∣74h3,and where fs∣7(6) and fs∣7(4) are the central difference representations of the 6th and 4th order derivatives of exact solution f(x) supported on 7 nodes of 2 adjacent FEs, respectively.The discontinuity of displacement derivativeu′j∣j+1(x) calculated above can be conveniently expressed in a nondimensional manner, which allows one to use it together with (10) in place of constant A. In such a case, it may be expressed as(16)p=1⟶Δ~sp=1≈-122f~s∣32,p=2⟶Δ~sp=2≈+144f~s∣54,p=3⟶Δ~sp=3≈-54166f~s∣76+32164f~s∣74.It should be realised that for higher degrees of approximation polynomialsp>3 or for nonuniform distribution of nodes within Fes, corresponding relations for the discontinuity of displacement derivatives u′j∣j+1(x) 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 p=1,2,3.The results on analytical calculations presented in Figure4 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 q that can be observed is directly related to the degree of approximation polynomials p and equal to q=p-1. So there are no frequency band gaps (q=0) in the calculated frequency spectra for the FEM numerical model based on the first degree of approximation polynomials p=1, while in the case of the third degree of approximation polynomials p=3 there are two frequency band gaps (q=2).Figure 4
Frequency band gaps of a 1D periodic medium, calculated analytically based on (10) and (16), for various degrees of approximation polynomials p=1,2,3 (c=5km/s, L=2m).The widths of frequency band gaps are also directly related to the degree of approximation polynomialsp. They appear for wavelengths λ equal to the FE length a or its multiples/fractions, that is, for p=2 at λ/a=2 or for p=3 at λ/a=1,2. However, it should also be noted that for multiple frequency band gaps associated with degrees of approximation polynomials p>2 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 FG1p=3 in comparison with FG2p=3.The results presented in Figure5 correspond to numerical calculations by the FEM for the same three degrees of approximation polynomials p=1,2,3. 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 f(x). 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 polynomialsp=1,2,3—a numerical model of 300 DOFs according to the elementary theory of rods [36], based on equidistant node distribution (c=5km/s, L=2m).It turns out that the behaviour observed in Figures4 and 5 stays closely correlated not only with the degree of approximation polynomials p, 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 FG1,…,FG4 appear at multiples of total number N 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 Figure6. It is assumed that the rod is made out of aluminium alloy of elastic modulus E=67.5GPa, material density ρ=2700kg/m3, and remains free—its both ends remain unsupported. The length of the rod is L=5m, while its height and width are the same a=b=10mm. For the assumed rod theory, phase cp and group cg velocity of propagating elastic waves are the same. For the given material properties, they are equal to c=5km/s. In such a case, longitudinal waves propagating within the rod require 1 ms to travel over length L of the rod from point P1(x=0) to point P2(x=L).Figure 6
Geometry of an isotropic rod of free ends.At first, the analysis of rod natural frequency spectra for various degrees of approximation polynomialsp 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 p increases the number of frequency band gaps, denoted as FG1,…,FG6. Their number is directly related to the degree of approximation polynomials p and is equal to p-1. As a consequence, for linear approximation polynomials p=1 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 polynomialsp—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 (c=5km/s, L=5m).Based on the results presented in Figure7, it may seem reasonable to reduce the degree of approximation polynomials p down to p=1 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 p. Additionally, it should be noticed that the degree of approximation polynomials p=1 corresponds to equidistant node distribution characteristic for the classical FEM.This behaviour is well illustrated by Figure8, in which relative errors ϵn of rod natural frequency spectra are presented for the same degrees of approximation polynomials p, which is the second part of the analysis. Relative error ϵn was calculated based on well-known analytical formula [34] for the type of boundary condition considered:(17)ϵn=fn-f^nf^n×100%,f^n=n2LEρ,n=1,…,DOF,where fn denote values of natural frequencies calculated numerically, while f^n is their corresponding analytical values.Figure 8
Relative errors of frequency spectra of an isotropic rod of free ends for various degrees of approximation polynomialsp—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 (c=5km/s, L=5m).It is well seen from Figure8 that the use of higher degrees of approximation polynomials results in a significant increase in relative error ϵn 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 p are characterised by much lower values of relative error ϵn.This error depends not only on the degrees of approximation polynomialp, 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 ϵ1 and ϵ2 were defined in a very simple manner as(18)ϵ^1=1N1∑n=1N1ϵn,ϵ^2=1N2∑n=1N2ϵn,2N1=N2=DOF.Table 1
Average errors of frequency spectra of an isotropic rod of free ends for various types and degrees of approximation polynomialsp—a numerical model of 736 DOFs according to the elementary theory of rods [36] used for calculations by the TD-FEM (c=5 km/s, L=5 m).
Polynomial type
Polynomial degreep
ϵ
^
1 [%]
ϵ
^
2 [%]
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.85The results collected in Table1 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 p errors ϵ^1 and ϵ^2 are smaller for approximation polynomials based on Legendre node distribution. As before, it should be emphasised that the degree of approximation polynomials p=1 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ϵ^1 decreases rapidly with an increase in the degree of approximation polynomial p. In the case of approximation polynomials based on Chebyshev node distribution, error ϵ^1 decreases 73 times, while in the case of Legendre polynomials it decreases 192 times. Lower errors ϵ^1 and ϵ^2 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 2p-2.Gauss quadrature employed in the case of Chebyshev node distribution is exact up to the degrees of approximation polynomials equal to2p 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 p-1 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 ϵn 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 (c=5km/s, L=5m).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 (c=5km/s, L=5m).


## 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 Figures7 and 9, force vibration responses of the rod were calculated assuming an excitation force F(f) variable within a frequency band from 0 Hz to 350 kHz and acting at point P1  (x=0), where also the amplitude q^ of rod vibration was assessed. The amplitude of the excitation F^ was 1 N. A small value of material damping was assumed, expressed by the loss factor η=10-4, to bound vibration responses of the rod nearby its resonant frequencies.Within the assumed frequency band frequency band gap FG3, 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 (c=5km/s, L=5m).It is evident that the presence of frequency band gap FG3 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 (c=5km/s, L=5m).During the analysis of time responses of the rod, various carrier frequenciesfc of the excitation force F(f) acting at point P1 (x=0) were examined with their spectra gradually closing to frequency gap FG3 [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 Figure13. In this figure, times t1 and t2 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 t1 and t2 are dependent on the excitation carrier frequency fc and their values are shown in Table 2.Table 2
Spans of propagating signal time windows for various excitation carrier frequenciesfc calculated analytically according to the elementary theory of rods [36] (c=5 km/s, L=5 m).
Carrier frequencyfc [kHz]
t
1 [ms]
t
2 [ms]
75
2.0
2.160
100
2.0
2.120
125
2.0
2.096
150
2.0
2.080Figure 13
Time responses of an isotropic rod of free ends for various excitation carrier frequenciesfc 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 (c=5km/s, L=5m).The results obtained clearly indicate that for excitation carrier frequencies distance from frequency band gap FG3, 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 FG3, these patterns become gradually deformed. In the case of the exaction carrier frequency fc 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 FG3 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 Figure15. 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 frequenciesfc 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 (c=5km/s, L=5m).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 Figures16 and 17 that present frequency spectra of the same rod obtained for the two-mode Mindlin-Herrmann theory of rods. The symbols S0 and S1 denote the first and second symmetrical (longitudinal) wave propagation modes, while cl and cs 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 (c=5km/s, L=5m).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 (c=5km/s, L=5m).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.Figures16 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 f^1, 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 f^1.Figure 18
Dispersion curves for the group velocitycg for the two-mode Mindlin-Herrmann theory of rod (cl=6086.1km/s, cs=3065.7km/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 FG4 just below the first cut-off frequency f^1 is much bigger than frequency band gap FG5. 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 f^1 is very clear.Finally, two important issues related to FE discrete numerical models should be addressed. The results of numerical calculations presented in Figure19 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 p employed. Moreover, the frequency spectrum is effectively divided into p-1 segments with p-1 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 (c=5km/s, L=5m).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 Figure20.Figure 20
Relative errors of frequency spectra of an isotropic rod of free ends for various dispersions of element lengthδl—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 (c=5km/s, L=5m).It was assumed that the lengths of FEs used for modelling the rod under consideration could vary within a certain limit, denoted asδl, but still these lengths are concentrated around a common average value l0=L/N. It can be noted from Figure 20 that an increase in the dispersion of rod element lengths δl slightly affects the observed error ϵn of rod natural frequencies. The strongest influence is observed in the case of the dispersion of rod element lengths δl equal to ±15%. In this case, the frequency band gaps FG1, FG2, and FG3 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 FG4 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.



---

*Source: 9651430-2016-01-05.xml*