Comparison of the Natural Vibration Frequencies of Timoshenko and Bernoulli Periodic Beams

This paper deals with the linear natural vibrations analysis of beams where the geometric and material properties vary periodically along the beam axis. In contrast with homogeneous prismatic beams, the frequency spectra of such beams are irregular as there exist enlarged intervals between some adjacent frequencies. Presented here are two averaged models of beams based on the tolerance modelling approach. The assumptions of classical Euler–Bernoulli and Timoshenko–Ehrenfest beam theories are adopted as the foundations. The resulting mathematical models are systems of differential equations with constant, weight-averaged coefficients. This makes it possible to apply any classical method of solution suitable for homogeneous beams, such as Galerkin orthogonalization. Here, emphasis is placed on the comparison of natural frequencies neighbouring the frequency band-gaps that are obtained from these two theories. Two basic cases of material and geometric property distribution in a periodicity cell are studied, and the natural frequencies and mode shapes are obtained for a simply supported beam. The results are supported by a comparison with the finite element method and partially exact solutions.


Introduction
Here, the considered periodic beams consist of repeated identical segments made from two or more different materials or segments made of a single material but with a variable cross-section dimensions or shape. Such structural elements possess certain distinguishing dynamic features, such as the existence of gaps in the frequency spectrum. The very mechanism of the gap occurrence in finite beams is clear. It is caused by the periodic distribution of the stiffness and mass properties that tends to favour certain vibration shapes and, as a result, the corresponding frequencies.
More precisely, as every mode shape consists of a finite number of half-waves, it is important if the mode shape nodes or antinodes are close or overlap with maximum or minimum stiffness and mass points. As a consequence, the frequency gaps for a simply supported beam are expected near frequencies for which the corresponding number of half-waves n is an even multiple or divisor of number of cells L/l (for an even number of cells) or L/l ± 1 (for an odd number of cells).
This paper focuses on linear-elastic beams with a periodic structure considered in the framework of the Euler-Bernoulli (abbreviated: Bernoulli) and Timoshenko-Ehrenfest (abbreviated: Timoshenko) theories. Generally speaking, the former of these theories has its application limited to analysis of slender beams, as it is a special case of the latter. The second theory takes into account both shear deformations and rotatory inertia effects, which makes it suitable for the analysis of moderately thick beams.
Many interesting remarks about these theories and their applications in the analysis of stepped beams can be found in the works [1][2][3][4][5][6], the first of which draws attention to the authorship of the so-called Timoshenko theory. There are many engineering fields in which multiple-stepped beams can be applied. Structures, such as aircraft wings, long span truss and beam bridges, suspended pipelines and rotatory unbalanced shafts can be approximately represented as such. For example, high buildings made from repetitive floors can be modelled as clamped beams with effective cross-section and equidistant lumped masses. There are other problems, such as periodically damaged structures, e.g., cracked reinforced concrete beams. In many of these domains, high-frequency vibrations are of crucial importance. Future applications are in the fields of vibration control and acoustic insulation.
Thus, it is desirable to develop efficient methods to make the analysis of low and high frequencies of these structures possible. The question arises: can an averaged model yield correct results for such non-continuous structures? The answer to this question can be pursued in various ways, including through numerical experiments as presented in this paper.
Dynamical problems of beams with periodic distributions of material and geometric properties are governed by differential equations with highly oscillating, often noncontinuous coefficients. There are various ways to tackle this problem. One group of these methods is based on discretization by means of the Finite Element Method, Finite Difference Method, or others. The second group strives to obtain simplified models, often based on homogenization. The resulting models replace the original, non-homogeneous structure with an equivalent homogeneous one with effective constant or slowly varying properties.
Among them, the theory of asymptotic homogenization [7] is one of the most popularized due to its mathematical rigorousness. However, these models often neglect the effect of the periodicity cell size, which is sufficient in the static analysis of such structures. In this paper, the tolerance-averaging technique (TA) [8] is applied in order to obtain averaged, non-asymptotic models of periodic beams. The resulting differential equations have constant, weight-averaged coefficients, some of which are dependent on the periodicity cell dimension.
In [9], Timoshenko's composite beam functions and Mindlin-Reissner thick plate theory were applied to a non-linear FE analysis of thin to moderately thick plates and slabs made of reinforced concrete. Non-classical, microstructure dependent Timoshenko beam models were analysed in [10][11][12], respectively. The effects of the moving load, boundary condition and material gradation of a bi-directional functionally graded beam on free and forced vibrations were examined in the framework of thin and thick beams in [13].
Periodically supported Timoshenko beams were analysed in [14]. Dynamic problems of periodic Timoshenko beams resting on a two-parameter elastic foundation were investigated in [15] by means of the weak-form quadrature element method. In [16], the problem of wave propagation in periodic Timoshenko beams on elastic foundations under moving loads, taking into account tensile and compressive axial load, was analysed.
Theoretical and experimental analysis of flexural vibration band gaps in Timoshenko beams with locally resonant structures was presented in [17]. In [18], the third order shear deformation theory was applied in the analysis of free vibrations of two directional functionally graded beams, taking into account various sets of boundary conditions. Nonlinear vibrations of Timoshenko beams using FEM were investigated in [19].
The cell-centre finite volume method was applied in the analysis of static and natural non-linear vibration analysis of functionally graded beams in [20]. Composites made of an isotropic elastic matrix that contain periodically placed inclusions or voids were analysed in [21].
A method of obtaining specified or extreme effective stiffness via topology optimization of the microstructure design was developed in [22]. Vibrations of composite structures, according to various beam theories, were considered in several papers. The beam can be modelled as a Timoshenko beam [23,24], Rayleigh beam [25,26] or Euler-Bernoulli beam [27].
Dynamic analysis can be extended based on the Winkler model [28,29], Pasternak model, a combination of both (the Winkler-Pasternak model), a nonlinear elastic model and fractional order viscoelastic model [30].
Researchers have presented methods to analyse structures resting on elastic foundations. The stability of periodic shells on elastic foundations subjected to external loading was considered in [31]. Papers [31,32] investigated the periodic behaviour of functionally graded plates and shells, respectively. Thin walled beams were analysed in [33], whereas the stress distributions and capabilities through a simple numerical example were demonstrated.
A theoretical study on the propagation of the flexural wave in the periodic beam on elastic foundation was presented in [34]. Moreover, the waves propagation was theoretically and experimentally investigated in straight beams with a periodic structure [35].
The aforementioned method was also applied to thermal problems [46][47][48] of periodic laminates, longitudinally graded materials, micro-heterogeneous hollow cylinders and cylindrical composite conductors, respectively. This work is a follow-up of earlier papers published by the author, together with other participants, on vibrations of periodic beams [49,50].
Although the results provided in these papers were satisfactory to some extent, the theory has been improved since then, mainly due to the use of a new class of weakly slowly varying functions in the analysis of shear-deformable beams. This contributed to an increase in the scope of application of the proposed models. An attempt was also made to compare the results of the proposed model with the exact results.
As stated before, vibration analysis of periodic beam type structures has recently received a remarkable amount of attention due to the importance and various applications of the subject. The tolerance-averaging approach is an effective analytical-numerical technique for problems that concern periodic structures. The main aim of this paper is to compare natural vibration frequencies and mode shapes of Bernoulli and Timoshenko beams using this approach.
The paper is outlined as follows. The theoretical background of the Bernoulli and Timoshenko beams is contained in Section 2. The tolerance approach and its fundamental assumptions are introduced and briefly discussed in Section 3. A brief derivation of the model equations in terms of both theories is presented in Section 4. In Section 5, the considered calculational problem is stated, and the solution methodology is described. In Section 6, the results and discussion are provided. The paper ends with our general conclusions.

Formulation of the Problem
We considered a beam made of linearly elastic material. The geometry of the beam is described in an orthogonal Cartesian coordinate system Oxyz. The Ox axis coincides with the axis of the beam, Figure 1. The inertial forces act in the direction of the axis Oz, and the Oxz plane of the load is, at the same time, the symmetry plane of the beam cross-section. With these assumptions, the problem becomes one-dimensional.
The region occupied by the beam is defined as Ω ≡ [0, L], where L stands for the beam length. We assume that the beam consists of many repetitive small elements, called periodicity cells. The periodicity cell is defined as ∆ ≡ [−l/2, l/2], where l L is the length of the cell. It is assumed that the beam thickness is of the order of the inhomogeneity period and that the deflections are small compared to the beam thickness.

Governing Equations of Euler-Bernoulli Beam Theory
The Euler-Bernoulli beam theory assumes that the sections perpendicular to an undeformed beam axis remain planar and normal to the axis after deformation. This theory also neglects the effect of rotational inertia. Let w = w(x, t) be the transverse deflection, EJ = E(x)J(x) be flexural stiffness, µ = ρ(x)A(x) be mass per unit length. Let ∂ k = ∂ k /∂x k be the k-th derivative of a function with respect to the x-coordinate; overdot stands for the derivative with respect to time.
The strain-displacement relations are assumed as follows: where κ stands for the beam axis curvature. Then, the relations between stress resultants (internal forces) and displacements can be introduced : where M is the bending moment. The strain and kinetic energy density per unit length of the beam are given by: The equations of motion can be obtained from the principle of stationary action A, formulated as: cf. [8], where the Lagrangian L = L x, t, w, ∂ 2 w,ẇ is given by the formula: and since the variations of unknown functions have to vanish at t = t 0 , t 1 , Equation (4) can be rewritten as: The differential equation of motion resulting from the above assumptions can be written as:

Governing Equations of Timoshenko Beam Theory
The Timoshenko beam theory includes both shear deformation and rotational inertia effects. This makes the theory applicable in the analysis of thick beams, sandwich composite beams, or beams subjected to high-frequency excitation when the wavelength is on the order of the beam thickness. In addition to the previous section, we introduce some new denotations: let θ = θ(x, t) be the average cross-section rotation and γ = γ(x, t) be the shear angle. Let kGA = k(x)G(x)A(x) be the shear stiffness and ϑ = ρ(x)J(x) be the rotational moment of inertia per unit length.
The strain-displacement relations are assumed as follows: as the cross-section rotation is no longer assumed to be equal to the deflection slope. Then, the relations between stress resultants (internal forces) and displacements can be introduced: where M, and Q are the bending moment and shear force, respectively. The strain and kinetic energy density per unit length of the beam can then written as: Again, equations of motion can be obtained from the principle of stationary action A (4); however, now, the Lagrangian (5) is a function of different arguments, L = L(x, t, w, θ, ∂w, ∂θ,ẇ,θ), and we obtain: The system of coupled differential equations of motion resulting from the above assumptions can be written as: The coefficients of Equations (6) and (12) are highly oscillating, often non-continuous functions of the x-coordinate. Subsequently, we present two averaged models of beams in order to circumvent this drawback.

Introductory Concepts and Basic Assumptions of the Tolerance Modelling
The tolerance modelling (or tolerance-averaging) technique (TA), cf. [8] is based on a set of intuitive heuristic concepts, such as tolerance relations, slowly varying functions and fluctuation shape functions. The most important of them are described here.
The operation of averaging over a region of periodicity cell is defined as: where y are coordinates in the local system associated with the cell, Ω ∆ is the region containing cell centres.
Let Ω be a regular region in R m and λ be a positive real number. Points x = (x 1, . . . , x m ) and y = (y 1, . . . , y m ) that belong to Ω are in tolerance determined by λ if the distance between those points is equal or less than λ. Now, let δ be a positive number. Real numbers µ, ν are said to be in tolerance determined by δ if the absolute value of the difference between these numbers does not exceed δ. These relations can be written as: A function F(·) will be called slowly varying, F ∈ SV R δ (Ω, ∆), of the R-th kind with respect to cell ∆ and the set of tolerance parameters δ = (λ, δ 0 , δ 1 , . . . , δ m ) if and only if the following conditions are satisfied: A function F(·) will be referred to as weakly slowly varying, F ∈ WSV R δ (Ω, ∆), if only the first of the above conditions is satisfied. The tolerance parameter λ is known a priori as a diagonal of the periodicity cell, and the δ k parameters can be determined a posteriori. Let h(·) be a λ-periodic highly oscillating function defined in Ω and on its boundary, continuous together with its gradients ∂ k h, k = 1, . . . , R − 1, and let it have piecewise continuous bounded gradient ∂ R h. The function is called the fluctuation shape function of the R-th kind h ∈ FS R δ (Ω, ∆), if it depends on λ as a parameter and satisfies the following conditions: for µ(·) being a certain positive valued λ-periodic function defined in Ω.
The tolerance-averaging approximation will be introduced only for the cases applicable in this paper. Let b(·), c(·) be arbitrary integrable λ-periodic functions (which usually are the medium physical properties) defined in Ω, and let us introduce the functions: For example, for C = 1, the tolerance approximations have the form: and the terms O(λ) are negligible. The micro-macro decomposition is based on the observation than the response of a periodic structure is periodic-like (periodic with respect to cell ∆ and tolerance parameters). Let us then decompose the transverse deflection and cross section rotation angle into their slowly varying and tolerance periodic parts (summation convention holds for superscript indices): where the decomposition of w(x, t) is valid for Bernoulli beam theory only, and decompositions of both w(x, t) and θ(x, t) are valid for Timoshenko beam theory. The new unknowns-the averaged transverse deflection, cross-section rotation and their fluctuation amplitudes-are weakly slowly varying functions of the C kind, respectively: and the corresponding fluctuation shape functions (FSs) are λ-periodic highly oscillating functions: which should approximate the deviation of the displacements in a cell from the average motion caused by the periodic structure. In the above relations, values C = 1, 2 define the classes of the unknown functions that are determined for the Timoshenko and Bernoulli beam theories, respectively.

Averaged Models
After substitution the micro-macro decompositions (19) into the Lagrangian (5), the next step of modelling is averaging (13) over an arbitrary periodic cell with approximations (18). The variations of the unknown functions are as follows: Now, the variation of the averaged action functional has the following form:

Timoshenko Beam Tolerance Model
Let us now introduce tolerance-averaged bending moments and shear forces for the Timoshenko beam model: The 2(1 + N) differential equations of the tolerance model for the macrodisplacements W(·), Θ(·) and their fluctuation amplitudes V A (·) and Z R (·) can be now written as: Finally, a system of differential equations with constant coefficients is obtained. The number of these equations depends on the number of introduced fluctuation shape functions. Additionally, some of the coefficients depend on the size l of the periodicity cell. From the principle of stationary action, we can also conclude the natural boundary conditions: Together with the averaged equation of motion, the following natural boundary conditions (for x = 0, L) with averaged coefficients are obtained:

Bernoulli Beam Tolerance Model
Let us now introduce tolerance weight-averaged bending moments for the Bernoulli beam model: This leads to a system of N+1 differential equations for macro-deflection and its fluctuation amplitudes: and natural boundary conditions: where W(x, t), V A (x, t) are the new kinematic unknowns. Together with the averaged equations of motion, the following natural boundary conditions (for x = 0, L) with averaged coefficients are obtained: It is worth mentioning that expressions (28) and (32) reduce to classic natural boundary conditions for a homogeneous beam.

Problem Statement
We consider an elastic beam of length L with l-periodic variation of geometric and material properties. These two types of properties variation are: variable height δ(x + l) = δ(x) and variable material properties E(x + l) = E(x), ρ(x + l) = ρ(x). The first case is related to a piecewise change of beam cross-section (stepped beam) for a beam made of a homogeneous material. Then, for steel the Young modulus E s = 205 GPa, the mass density ρ s = 7850 kg/m 3 , the Poisson's ratio ν s = 3/10 are assumed. The cross section is rectangular In the second case, a beam consists of repetitive sections of aluminium and steel; therefore, different Young moduli and mass densities are considered (bi-material beam). For aluminium, the Young modulus E a = 69 GPa, the mass density ρ a = 2700 kg/m 3 , the Poisson's ratio ν a = 3/10 are assumed. The cross section is rectangular: The length of the beam for both cases, is L = 1.0 m. The shear coefficient was assumed equal to k = (10 + 10ν)/(12 + 11ν).
A fragment of the beam and a single periodicity cell are illustrated in Figure 2. In all of the cases, the length of the cell length was assumed to be equal; l = L/10. The analysis was performed for seven values of the saturation parameter α = 1/8, 1/4, . . . , 3/4, 7/8, which stands for the relative length of the central segment of the cell. l l l

Effective Properties
Determining of the fluctuation shape functions is a crucial point in calculation of the averaged coefficients of the Equations (26) and (30). These functions have to, at least approximately, represent the eigenmodes of a periodicity cell. From previous works, we concluded that the best results were obtained when analysing a combination of two cells. A finite element model of a two-cell assembly with periodic boundary conditions was used to obtain the eigenmodes. The minimal number of elements per unit cell was set by a prior convergence study. The first four of the fluctuation shape functions h A for stepped and bi-material beam are displayed in Figure 3.

Natural Vibrations Analysis
In order to obtain a system of algebraic frequency equations, the Galerkin method is applied. The trial functions are assumed in the form of truncated trigonometric series. In the case of a Timoshenko beam, these are: where the x-dependent functions X m and Y m have to satisfy the appropriate boundary conditions. The macro-displacements solutions for a simply supported beam were assumed as X m (x) = sin(mπx/L) and Y m (x) = cos(mπx/L). For the fluctuation amplitudes V A and Z A , which come in pairs, one has to carefully examine the boundary conditions. The fluctuation shape functions for a symmetric unit cell are of two kinds. If the fluctuation shape function h A is symmetric, its corresponding p A function is antisymmetric with respect to the centre of the cell, and vice versa. In that case, X V A m (x) = X Z A m (x) = sin(mπx/L); and in the opposite case, X V A m (x) = X Z A m (x) = cos [(m − 1)πx/L]. In case of the Euler-Bernoulli theory, the rotation of cross section is dependent of the deflection, and naturally unknowns Θ(x, t) and Z R (x, t) drop out of the equations.
Inserting the assumed solutions to either the Timoshenko (26)

Exact Models and Their Numerical Solutions
Let us divide the beam into n sections wit constant geometrical and material properties, so that EJ = const, kGA = const, µ = const, ϑ = const. Then, in every section, the equations of motion according to the Timoshenko theory have the form: Now, we introduce the nondimensional coordinate ξ i = x i l i for i-th segment, i = 1, . . . , n. The relation between derivatives can be written as: Solutions to the obtained differential equations can be sought in the form: where the coefficients are as follows: The internal forces in i-th section of the beam are equal to: Now, one has to define the continuity conditions of displacements and internal forces on the interfaces between the neighbouring sections: Moreover, the boundary conditions have to be defined on the beam ends. For a simply supported beam these conditions are: From the above relations, we obtain a set of homogeneous algebraic equations for the constants C 11 − C n4 : The condition for existence of nontrivial solutions is as follows: det(A(ω)) = 0, (44) from which we obtain a transcendental equation, which was solved in two steps. First, the roots were isolated by linear search. Next, the secant method was applied to determine the roots with precision The solution procedure has to be simplified in order to analyse this problem in the framework of Bernoulli theory. Assuming that the shear angle γ is now equal to zero, the average cross-section rotation θ = ∂w, we obtain Equation (7), bearing in mind that now the bending stiffness is constant. The solution for the deflection w given by first of the expressions in Equation (36) is still valid. Now, in the expressions (38) and (39) i . The coefficients in Equation (37) are now equal, a 2 i = b 2 i = λ 2 i . The boundary and continuity conditions remain the same, but the bending moment and shear force in Equation (40) have to be written as: and the general form of Equations (43) and (44) remains unchanged.

Results and Discussion
In this section, a numerical study of natural frequencies of considered Bernoulli and Timoshenko beams will be carried out. The influence of beam stiffness and mass distribution will be investigated in relation to the saturation parameter α. First, low order frequencies ω 1 -ω 5 for the two considered cases, for various saturation parameter α = 1/2, according to the proposed and the exact models, are investigated.
The results are displayed in Table 1. The discrepancies between the averaged and exact models increase slowly with increasing mode number, and the agreement is remarkable. The relative differences are less than 0.0075%. Then, frequencies ω 1 -ω 5 for the two considered cases, for various values of saturation parameter α, according to both models, are investigated. The results are displayed in Tables 2 and 3 for a stepped beam according to the Timoshenko and Bernoulli theories, respectively. The results for a bi-material beam are shown in Tables 4 and 5. The results of the TA and FE models are in excellent agreement, as the relative difference between these models does not exceed 0.02%.
It can be seen that the low frequencies increase with increasing α parameter value for a stepped beam. In the case of a bi-material beam, these frequencies are lowest for α = 1/2 and increase when this parameter decreases towards its minimum value and when it increases towards its maximum value. For these low frequencies, the relative difference between Bernoulli and Timoshenko model is less than 5% for the stepped beam and less than 3% for the bi-material beam.   The frequency spectrum for a stepped beam is presented in Figure 4 for Bernoulli and Timoshenko beam models for the saturation parameter α = 1/2. The first frequency gap is placed between the 9th and 10th natural frequency and is equal to ∆ω =14,493 rad/s according to the Bernoulli beam model and ∆ω =10,588 rad/s according to the Timoshenko beam model. The second frequency gap can be found between the 20th and 21st natural frequency and is equal to ∆ω =34,728 rad/s for the Bernoulli model and ∆ω =17,594 rad/s for the Timoshenko model.  Similar results are presented in Figure 5 for a bi-material beam and for the saturation parameter value α = 1/2. There is no distinct frequency gap near the 10th frequency for either of the models. The difference between the 20th and 21st frequency is equal to ∆ω =34,378 rad/s and ∆ω =22,418 rad/s for the Bernoulli and Timoshenko model, respectively.   For the stepped beam case, the maximum difference occurs between the 9th and 10th frequency for α = 3/4 for both models, Figure 6. For higher frequencies, the maximum difference occurs between the 19th and 20th frequency for α = 7/8 for the Bernoulli model; however, according to the Timoshenko model, it is between the 20th and 21st frequency for α = 1/2, cf. Figure 7.
For the bi-material beam case in the lower frequency spectrum, the largest gap occurs between 10th and 11th frequency for α = 1/4 for both models, Figure 8, and there is no distinct gap near α = 1/4, which confirms the data in Figure 6. For higher frequencies, the maximum difference occurs between the 20th and 21st for α = 1/2 according to the Bernoulli model, and for α = 5/8 according to the Timoshenko model, cf. Figure 9.
The general conclusion that can been formulated from the considerations is that the frequencies obtained from the Timoshenko beam model are lower then the calculated from the Bernoulli model, and the differences increase with increasing the number of eigenmode half-waves.
Considering the differences between the results obtained from the compared theories, let us introduce the relative differences r n = (ω (n) T , corresponding to the frequency number n.
The growth of these differences, as a function of the frequency number, is approximately quadratic; at the same time, it is nearly linear as a function of frequency. This can be said about both the stepped and bi-material beam cases. For example, the relative difference for the stepped beam varies from r 1 = 0.1% for α = 1/8 and r 1 = 0.2% for α = 7/8 to r 21 = 39.3% for α = 1/8 and r 21 = 67.7% for α = 7/8. In the case of a bi+material beam, these values vary from r 1 = 0.1% for all values of α to r 19 = 31.9% for α = 1/8 and r 19 = 45.8% for α = 7/8. The situation is therefore, unsurprisingly, similar to that of uniform beams of constant height.
On the other hand, it is expected that the heterogeneity at the cell level should have some impact on the results. With fixed proportions of cell sections' stiffness, the parameter α should be the crucial parameter. Let us first analyse the results for the stepped beam.
When increasing this parameter from its lowest to the highest value, the relative differences increase as well. However, up to the fifth form of vibration, these differences do not exceed 5%. For the sixth form, these values range from r 6 = 3.8% for α = 1/8 to r 6 = 6.9% for α = 7/8. For the 10th mode, the differences range from 11% to 21%. Generally, for the higher frequencies, the relative frequencies increase faster.
For the bi-material beam, the distribution of differences is of a slightly different nature. Initially, as in the previous case, these differences are less than 5% for the first five modes of vibration. Up to the eighth frequency, the discrepancies do not exceed 8% and are basically independent of the parameter α. For higher modes, the relative differences tend to grow from their minimal values for extreme values of α to their maximal values for intermediate values of this parameter. However, the changes in the discrepancies are much more moderate than for the stepped beam. For the 11th mode, they are equal to r 11 = 14.2% for α = 1/8 through r 11 = 16.7% for α = 5/8 and to r 11 = 14.7% for α = 7/8.
In Figures 10 and 11 the 9th, 10th and 11th mode shapes of the considered beams are displayed. The shapes are shown for the left half of the beam, x ∈ (0, L/2) for the stepped beam ( Figure 10) and for the bi-material beam ( Figure 11). It can be seen that the TA and FE results are in good agreement, but the deflections obtained from the Bernoulli beam model differ slightly for some cases, especially near the beam supports, as the shear forces are maximal there. (a) (b) (c) Figure 11. The 9th (a), 10th (b) and 11th (c) natural vibration modes presented for the left half of the beam, for a bi-material beam for the TA and FE Timoshenko beam models and for the TA Bernoulli model.

Conclusions
In this paper, the basic dynamic properties of beams with periodically varying geometric and material properties were investigated. The periodic variation of geometric and material properties produces the band gaps in the frequency domain. It was shown by comparison with FE results that the proposed tolerance-averaging technique for differential operators leads to reliable homogenized models of beams.
As some of the averaged coefficients of obtained differential equations depend explicitly on the unit cell length, the positions of gaps in the frequency spectrum were determined properly. They can also be predicted by considering the order of the vibrational forms of the periodicity cell. We confirmed that, for the considered boundary conditions, the trial functions used in the Galerkin method solutions can be adopted from solutions of a homogeneous beam.
As the second result, we demonstrated, through comparison between both models, that taking the Euler-Bernoulli assumptions should be considered with caution. In the higher frequency range, the simplification of neglecting the shear susceptibility and rotational inertia terms leads to unacceptable discrepancies in relation to the Timoshenko model, which was proven by many researchers to be more accurate.
Applications of the so-called technical theories of beams have their limitations. The Bernoulli model's application possibilities are limited to the analysis of rather slender beams. More precisely, it is required that the half-wave length of any eigenmode is significantly greater than the cross-section height. The Timoshenko model, concerning shear deformability and rotatory inertia as well, has the potential to yield correct results for thick beams.
Let us be concerned only with the technical 1-D theories of continuous beams here. The main factors that matter are the half-wave length of the n-th vibration mode λ n = L/n and the length-to-height ratio of the beam η = L/h. It was shown in [51] that the possible scope of application of the Timoshenko model is rather significant.
Comparing the results obtained from this model and a 2-D FE model, the half-wave length can be of the order of the beam thickness. Thus, the number of eigenfrequencies and eigenmodes is limited by the ratio of the mode half-wave length to the beam thickness. It was assumed here that the safe range of applicability of Timoshenko's theory is for the ratio λ/h max greater than 1.25.
A shortcoming of this paper is the modelling assumption that the bonds between neighbouring sections of the beam are ideal. Especially in stepped beams, there are stressfree areas of the beam adjacent to the height jump. Therefore, bending stresses do not occur over the entire height of the deeper section. This could be taken into account by assuming a functional change in the cross-section height in the vicinity of these points of transition between the two sections. The proposed models allow for the introduction of such an assumption. However, a significant qualitative change in the results is not expected.

Informed Consent Statement: Not applicable
Data Availability Statement: The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The author declares no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: TA Tolerance Averaging FE Finite Element