Nonlinear dynamics of a reduced multimodal Timoshenko beam subjected to thermal and mechanical loadings

Large amplitude vibrations of a Timoshenko beam under an influence of temperature are analysed in this paper. In the considered model the temperature increases instantly and the heat is uniformly distributed along the beams length and cross-section. The mathematical model, represented by partial differential equations takes into account thermal and mechanical loadings. Next, the problem is reduced by means of the Galerkin method, considering the first three natural vibration modes of a simply supported beam in the both ends. The influence of the temperature on amplitudes and localisation of the resonance zones and stability of the solutions is studied numerically and analytically by the multiple time scale method. The bifurcation points, existence of unstable lobes and transition from regular to chaotic oscillations are shown.


Introduction
The beams are fundamental structural elements with application in many branches of the industry. Frequently, these structures are subjected to dynamic loading leading to large amplitude vibrations. Large vibrations introduce a geometrical type of non-linearity that influences the dynamic behaviour of a structure. In this case the structures stiffness, and consequently the resonance frequencies and mode shapes, are amplitude dependent. Linear and nonlinear vibrations of beams have been deeply investigated for many years and were reviewed, for example, in the books of Nayfeh and Mook [1] and of Nayfeh and Pai [2]. Nonlinear vibrations of the Euler-Bernoulli or shear deformable beam models have been studied there and the influence of the nonlinear terms on the bifurcation scenario and possible resonances have been discussed in details. The advanced composite beam theory has been presented in [3]. The beam models considered various configurations of lamina with reinforced fibers orientation, closed or open cross-section shapes. The beams are also used to model rotating blades dynamics, for example blades of a helicopter rotor [4,5] or blades of turbines. Many authors use the classical FEM [6,7] and semi-analytical methods [8,9] to study this problem.
In most of the analysis the environmental conditions are neglected. One of very important factors which has to be considered is temperature, which may vary in high ranges in real mechanical or aerospace applications. Temperature variations can and do affect substantially the vibration response of a structure. Thermal loads introduce stresses due to thermal expansion, which lead to changes in the modal properties. The basic problems of the thermoelastic vibrations can be found in the books of Boley and Weiner [10], Nowacki [11] and Thorton [12].
Although the temperature and elastic behaviours are in fact coupled [12,13] for thin structures it is often reasonable to assume that the temperature distribution is independent of the deformation or that the structure gets the elevated temperature instantly. This approach is widely used to model the thermoelastic behaviour of structures. The geometrically nonlinear vibrations of structures at the elevated temperature are studied by many authors as [14][15][16], etc. In [17] and [18] thermomechanical, geometrically nonlinear vibrations of plates and beams, correspondingly, are studied. The authors found a very reach nonlinear dynamic behaviour of the system including, periodic, quasi-periodic and chaotic oscillations. A thermomechanical model of the vibration of a Timoshenko beam after its one mode reduction is studied by multiple time scale method in [19].
In the present work the study is extended by using three mode reduction of the beam's model for thermoelastic vibration. The goal of this paper is to show specific dynamic phenomena of geometrically nonlinear vibrations of a Timoshenko beam subjected to thermal and mechanical loadings. The phenomena, such as bifurcations, non-periodic or chaotic oscillations which arise due to varying temperature are taken into account in the study.

A model of a Timoshenko beam under thermo-mechanical loadings
The considered structure is a beam made of elastic composite material subjected to thermal and mechanical loadings. The beam orientation together with coordinates and indicated length l, thickness h, and width b, is presented in Fig. 1. The mathematical model of the Timoshenko beam presented in Fig. 1 has been derived in papers [13,16]. The dimensionless equations of motion have the form: where u, w, w are dimensionless displacement field coordinates, a = 12l 2 /h 2 , b = kG/E, k-shear correction factor, G, E-respectively shear and Young modulus, p-external mechanical loading. Also linear damping terms with damping coefficients d 1 , d 2 have been added to the model.
a T is a linear coefficient of a thermal expansion, and where h (x, z, t) is a space and time dependent dimensionless function of a temperature distribution.
Assuming that the distribution of temperature along x and z axis is constant: h (x, z ) = const., we get: where DT is a difference between the reference temperature and the current temperature. According to the paper [19] longitudinal displacement can be found from the first equation of the set (1) and the model can be simplified to two partial differential equations, having the dimensionless form: We note that the dimensionless displacement of the beam is expressed versus beam's length w = w * /l, where w * is the displacement in physical units. It has been assumed that mechanical loading p(x, t) is distributed along x axis and is expressed by a function of space and time.

Multimodal model reduction
The model of a beam has been reduced from partial differential equations (PDE) to ordinary differential equations (ODE) by means of modal projection and then by using Galerkin's orthogonalisation method. For this purpose the generalized displacements vector is expanded as a sum of the product of the quasi-normal modes w n (x), w n (x) and the time dependent functions q n (t): where N f is a number of assumed modes. Substituting Eq. (5) into (4) we get where According to the Galerkin's procedure the quasinormal modes should satisfy geometrical and dynamical (natural) boundary conditions. For the projection we take linear modes of beam natural vibrations reported in ''Appendix 1''. By using the fact that w n and w n are solutions of the eigenvalue problem we obtain: Àx 2 n w n q n ðtÞ À d 1 Multiplying (6) 1 by mode w m , and (6) 2 by mode w m (x ), then summing up both equations we have Then, integrating (9) over the beam length, invoking the orthogonality condition, and assuming proportional damping the equations are transformed into the form: where, x n is the nth natural frequency of the linear undamped Timoshenko beam and n n is a dimensionless modal damping coefficient. Let's consider the first three modes of the expansion (5), N f = 3. Applying the formulae (5)- (11) we get the equations of motion Coefficient K is calculated from (7) and has the form where The coefficients B n (n = 1, 2, 3) are defined in ''Appendix 1''. External loading has been assumed as distributed periodic force where P(x) is a space function of distributed loading and f(t) is a function of time varying loading which is accepted as a time periodic function with f a amplitude and X frequency of external loading. Substituting (13) and (15) into (12) we get a set of nonlinear equations The parameters p 1 , p 2 , p 3 for external loading distributed according to the normalised mode shape take definitions The values of coefficients of Eq. (16) are determined for data reported in ''Appendix 2'' for a case of a symmetric cross-ply laminated beam composed of 20 orthotropic layers.

Analytical solutions
Differential equations of motion (16) are nonlinearly coupled by cubic terms. Therefore, in order to solve the problem analytically the approximate method of multiple time scales is used [21]. The set of Eq. (16) is reordered by introducing a formal small parameter e, where, l 1 = 2n 1 x 1 , l 2 = 2n 2 x 2 , l 3 = 2n 3 x 3 . The formal small parameter e is used for grouping small terms on the right side of Eq. (18), thus the coefficients are defined as: The solution for q 1 (t ), q 2 (t ), q 3 (t ) is assumed in the form of a series of a small parameter e Generalised q i coordinate (i = 1, 2, 3) is expressed respectively in the zeroth, first and second order perturbations, q i,j (T 0 , T 1 , T 2 ), (j = 0, 1, 2). Dimensionless time is also expanded in a series of a small parameter where T 0 , T 1 , T 2 , are respectively the fast and slow time scales. The first and the second time derivatives are now defined with respect to the introduced time scales: The operator D m n ¼ o m oT n denotes the mth order partial derivative with respect to the nth time-scale.
Considering the three mode reduction we may expect the main resonances occurring around the natural frequencies x 1 , x 2 and x 3 . The analytical solutions are determined around these resonance zones.
The solutions are sought around the natural frequencies of the beam, thus excitation frequency X has to satisfy the condition where r i is the frequency detuning parameter of ith resonance zone (i = 1, 2, 3). The methodology of the analytical solution determination will be demonstrated for the first resonance zone around x 1 . For the other two cases the procedure is identical, just coordinates q 2 and q 3 play dominant role around the second and the third resonance, respectively.
In further analysis we exclude the internal resonance case. Therefore the ratios of the natural frequencies of the system are assumed to be incommensurable numbers. In order to excite a selected mode directly, we accept external excitation as: p 1 [ 0 and p 2 = 0, p 3 = 0. In such a case, near the first resonance zone the solution is manifested mainly by the first coordinate. The coupling with the other coordinates occurs in the higher perturbation orders. Therefore solution of Eq. (23) has the form is the imaginary unit, A 1 is the complex amplitude and " A 1 is complex conjugate. Substituting solutions (26) into (24) and, after grouping the terms in proper exponential functions, we get where cc means complex conjugate functions to those written in the right side of Eq. (27) and ST 1 represents secular generating terms. In order to avoid secular terms in the solution we require ST 1 to be zero The particular solution of (27) has the form Substituting solution (29) into (25) we get NST 1 , NST 2 are nonsecular generating terms and they are not directly reported here, and ST 2 are secular generating terms of the second order which should vanish Taking into account the solutions (26) and (29), expressing the complex amplitudes in the polar form using expansion (19), after transformation to the trigonometric form, we obtain solutions in the zeroth and the first order approximation Applying the so called reconstitution method [21], on the basis of (28) and (31) we formulate the modulation equations for the complex amplitudes A 1 as the first order ordinary differential equation Expressing complex amplitude A 1 in the polar form (32) and then separating the real and imaginary parts, we get the so called modulation equations for amplitude a 1 and phase / 1 Amplitude a 1 and phase / 1 can be found from the modulation Eq. (35), or for a steady state case from algebraic equations equalling the derivatives to zero, da 1 dt ¼ 0; d/ 1 dt ¼ 0: Then, we can find the resonance curve by determining sin/ 1 and cos / 1 and after some algebraic manipulations we get a polynomial type equation depending on the system's parameters: where z = a 1 2 . Coefficients b k , k = 0, … 7, are reported in ''Appendix 3''. On the basis of Eq. 36 we can study the influence of selected parameters on the beam's response. Having amplitude a 1 the approximate solution can be obtained from (33).
Stability of the solution is determined by analysis of the modulation Eq. (35) which can be written in the shorter form Perturbing Eq. (37), considering a linear part of the power expansion and then subtracting perturbed and unperturbed equations we get a set of first order linear differential equations in perturbations where d means perturbation (variation) of the variable. The stability depends on a real part of the eigenvalues of the Jacobi matrix The solution is unstable if at least one of the roots has real part positive. The solutions near the second and the third resonance zones are found following the same procedure.

Numerical results
Numerical calculations have been carried out for data presented in ''Appendix 2'' for a symmetric cross-ply laminated composed of 20 layers composite beam. Amplitude and frequency of external load are varied in order to demonstrate essential nonlinear phenomena around the resonance zones or bifurcation points.
The resonance curves obtained from analytical Eq. (35) are presented in Fig. 2. The curves are plotted for fixed temperature DT ¼ 20 and three selected levels of excitation p i = 10 -7 ,the so called modulation p i = 10 -6 , p i = 5 9 10 -6 , i = 1, 2, 3. The system is excited around the resonance assuming that the harmonic excitation corresponds to the excited mode, i.e. around the natural frequency x i only excitation p i is activated. For a small level of excitation amplitude p i = 10 -7 resonance curves are similar to a linear beam model. While increasing the amplitude the curves exhibit stiffening effect around three considered resonances. Assuming the same intensity of excitation, the resonance curves around the third natural frequency get the smallest amplitudes comparing to two other cases (see Fig. 2c, d). The shape of the lowest resonance curve presented in Fig. 2d is shown in Fig. 2e in an enlarge scale. Because the natural frequencies of the system are well separated there is no internal resonance in the structure, the coupling between modes exists only due to nonlinear geometrical terms and for small oscillations is not visible. The nonlinear dynamic phenomena arise mainly near the first resonance zone and for relatively large oscillations. The first observation is that there is a loss of stability at the beginning, on the left side of the resonance curve in Fig. 2a for p 1 = 5 9 10 -6 . The zoom of the beginning of this curve is visible in Fig. 2b.
The influence of the temperature can be determined from Eq.(35) (or 36). For the steady state we equalize these equations to zero and determine amplitude for selected values of temperature. To show the temperature influence we plot the resonance curves for three resonance regions and three selected temperatures values DT ¼ À70; DT ¼ 0 and DT ¼ þ70; We observe ( Fig. 3) that the increased temperature shifts the resonance region to lower frequencies (to the left side) and for the negative value the curve is moved to higher frequencies (to the right side), for all considered resonance regions. Apart from this the elevated temperature increases oscillation amplitudes. This phenomenon is visible around the first resonance zone in particular (Fig. 3a). Moreover, one can notice that Fig. 2 Resonance curves for different excitation level, p i = 10 -7 , p i = 10 -6 , p i = 5 9 10 -6 , i = 1, 2, 3, around a the first, c the second and d the third natural frequency; b zoom of the curve p 1 = 5 9 10 -6 , e zoom of the curve p 3 ¼ 10 À7 ; DT ¼ 20 unstable part of the curve is extended comparing to the temperature DT ¼ 20 in Fig. 2a.
The analytical results have a very good agreement comparing them with the ones obtained by direct numerical integration. In Fig. 4 a comparison of the analytical solution (solid line) with direct numerical integration of the original ODE (Eq. 16) by Runge-Kutta method is presented. We see a very good agreement for relatively large oscillations. A difference is observed near low frequency, on the left hand side of the curve around frequency X ¼ 0:01: Of course, the analytical method works well only if oscillations are not very large. For large oscillations we apply continuation method by Auto package [22], in order to discover bifurcation points and new nonlinear phenomena.
In the further analysis we focus around the first resonance zone considering large oscillations. We increase meaningfully amplitude of external excitation, p 1 = 4 9 10 -5 . In spite of this fact, the solutions obtained analytically (Fig. 5a) and numerically (Fig. 5b) does not differ essentially. The main difference is visible on the borders of the resonance. In contrast to the analytical solution the resonance curve obtained from direct numerical integration of ODEs is unstable, close to the maximal amplitudes. The zoom z 1 of this region is shown in Fig. 5c and more details of the zoom z 3 in Fig. 5d. We see that close to the top of the resonance curve, the solution loses stability and an additional unstable branch arises (Fig. 5d) which does not occur in analytical solutions. Moreover, on the left side about X ¼ 0:02 unstable solutions of large Fig. 3 Resonance curves for different temperature levels, DT ¼ À70; DT ¼ 0; DT ¼ þ70; around a the first p 1 = 5 9 10 -6 , p 2 = 0, p 3 = 0, b the second p 1 = 0, p 2 = 5 9 10 -6 , p 3 = 0 and c the third natural frequency p 1 = 0, p 2 = 0, p 3 = 5 9 10 -6 amplitude are obtained analytically which do not have confirmation by direct numerical simulations (see left side of Fig. 5a and zoom z 2 in Fig. 5f).
Comparing the one mode reduction given in [19] with the present study based on the three mode reduction we may conclude that for periodic excitation distributed according to the selected mode the involvement of the rest modes is small. For small vibration amplitudes the higher modes involvement can be neglected. A difference however, appears for larger oscillations. The course of the resonance curve for the one mode or three modes reduction is similar as presented in Fig. 5b but for the one mode reduction the top of the resonance curve (zoom z 1 ) is stable. This part of z 1 zone obtained for the one mode reduction is presented in Fig. 5e. Of course, we may expect that the involvement of a larger number of modes will be crucial if the structure is tuned to the internal resonance condition which is excluded from the present study. Also the multimodal model allows studying the system for a case of multimodal mechanical or thermal loadings.
We may notice that for low frequencies there are untypical lobes which are partially unstable (Fig. 5f). We may expect that this part of the curve is sensitive for the temperature change. The temperature influence is presented in Fig. 6. The lobes are almost stable if the temperature is low (DT ¼ À50). Elevated temperature (c) DT ¼ 0; (e) DT ¼ 50; (g) DT ¼ 70 makes the lobes unstable. Moreover, for temperature DT ¼ 70the lower branch on the right side becomes unstable too. Time histories (for one period) of the solution q 1 for temperature DT ¼ 70 and selected excitation frequency, starting from low frequency and then increasing it, is shown in Fig. 7. For the low frequency (Fig. 7a-c) the beam response includes more than one harmonic. Approaching the resonance zone the response is dominated mainly by a single harmonic corresponding to the excitation frequency ( Fig. 7d-g).
Increasing the temperature above a certain threshold two new branches occur in the resonance curve (Fig. 8). Close to the branch point these additional solutions are unstable, but they become stable while the frequency is increasing. It means that in a wide frequency domain there are five possible solutions, three of them are stable and two unstable. The solutions located on the left hand side of the curve are mainly unstable (see the zoom of low frequency lobs in Fig. 8b). The system looses stability in this region and we observe chaotic oscillations there (Fig.  9a). Depending on the initial conditions however, another regular attractor coexists together with the strange chaotic one. This is clearly visible in Fig. 9b in which different basins of attraction of regular (two points) and chaotic motion are present.
The elevated temperature changes dynamics of the considered beam from regular to chaotic. The chaotic attractor can be quenched by the decrease of the temperature. In the bifurcation diagram in Fig. 10 the reduction of the chaotic attractor is presented. Varying temperature from DT ¼ 100 till DT ¼ 70 we observe declining chaotic attractor with small periodic windows, and about DT ¼ 74 there is a transition to subharmonic regular motion of 9T-period and for about DT ¼ 70 there is a transition to subharmonic regular 2T-period motion. The scenario is presented in Poincaré maps in Fig. 11. There is a visible jump from chaotic to regular dynamics due to small temperature variation from DT ¼ 75 to DT ¼ 74:

Conclusions and final remarks
Dynamics of a Timoshenko beam model under an influence of the elevated temperature is analysed in this paper. The model described by partial differential Fig. 4 Resonance curve around the first natural frequency for excitation level, p 1 = 5 9 10 -6 ; solid line analytical solution obtained from ME, dots direct numerical integration of ODE  equations has been reduced to ordinary differential equations and the detailed study has been carried out for the first three vibration modes. On the basis of analytical solutions it has been shown that the system exhibits stiffening effect in the resonance curves for relatively large amplitude of excitation. Moreover, the instability of the solution, observed by unstable lobes, located on the left side of the resonance curve has been found. A particular attention in the study has been paid to the influence of the temperature on quantitative and qualitative changes in the beam response. It has been shown that elevated temperature shifts the resonance zone in the direction of lower frequencies and increases amplitude of oscillations. This phenomenon occurs around all three considered resonance areas, however this is the most exposed near the first natural frequency. An increase of the temperature makes unstable not only the lobes occurring on the left side of the resonance curve. Also the instability is observed on the right site of the curve for high frequencies. Above a certain temperature threshold two additional branches arise on the resonance curve. Furthermore, the elevated temperature for low excitation frequency may transit the beam into the chaotic oscillation region. Decreased   These characteristics correspond to a glass-epoxy composite material. The effective Young modulus of the beam is: E ef = 41.92 GPa.
The effective Young modulus have been computed considering layers having E 1 with orientation 0°and layers having E 2 with orientation 90°. Using the simple summation formulae: we obtain the equations of composite Timoshenko beam in the form: To obtain E ef in the usual form of the Timoshenko beam the coefficient D is divided by b