Transverse vibrations of a beam under an axial load: minimal model of a triangular frame

The paper is devoted to transverse in-plane vibrations of a beam which is a part of a symmetrical triangular frame. A mathematical model based on the Hamilton principle, formulated for large deflections of the beam subjected to dynamic axial excitation, is presented. An approximate nonlinear ordinary differential equation for the vibration amplitude is derived by means of the Galerkin method. Dynamics of the system is studied numerically for the two cases: harmonic and pulsating load. Various values of the model parameters, including the excitation amplitude and frequency, are considered. The amplitude is taken to be below or above the static critical load. The regions of stable and unstable solutions are determined in parameter planes by evaluation of the maximal Lyapunov exponent. The results are compared to the case of the abbreviated (linear) dynamical system. The stable and unstable beam responses are analysed and classified.

Mechanical system to be considered: a the full triangular frame, b the cross-beam isolated from the system

Mathematical model
Consider a symmetrical triangular frame (see Fig. 1a) whose arms (1) of length L 1 are connected via a crossbeam (2) of length L 2 . The arms have a rectangular cross section (b × c). The cross section of the beam, in turn, has the shape of a circular ring: its outer and inner diameters are denoted by d 1 and d 0 , respectively. The vertex C of the frame is fixed while the others, A and B, are simply supported out-of-plane and subjected to in-plane load.
Such a system can be regarded as a simple representation of a brake triangle that has been produced on a large scale for the railway industry. It is an important component of bogie brake systems suitable for various types of wagons. The triangle transmits the braking force from a brake cylinder to brake pads that, in turn, mate with rail wheels [15].
In paper by Magnucki and Milecki [14], one can find, among others, the relations between the internal forces in the cross-beam, i.e. the axial force N 2 and the bending moment M b2 , and the load F. Now, in-plane vibrations of the system are considered, i.e. the load is assumed to be a function of time, F = F(t). To simplify dynamic analysis, the beam is isolated from the structure as can be seen in Fig. 1b. The beam is treated as simply supported but subjected to the axial force N 2 (t) and opposite end moments M 2 (t). For simplicity energy dissipation is neglected.
Behaviour of the beam is described by the transverse deflection w(x, t). The axial displacement of a given point is u(x, z, t), where z denotes the distance from the neutral surface. A small element of the beam is illustrated in Fig. 2. In order to describe vibrations of the system, we take an approximate energy approach.
In case of large deflections, the normal strain can be written as For a linearly elastic material, the strain energy of the system is given by in which E is Young's modulus and are the area of the cross section and its moment of inertia about the neutral axis. It is assumed that the rotary inertia of the beam can be neglected, thus the kinetic energy has the form where ρ denotes the material density. Finally, the work done by the load is expressed as Note that the effect of the moments M 2 are not included at this stage. The application of Hamilton's principle leads to the following partial differential equation of the transverse beam motion: The first three terms of this equation are the same as in case of the linear model based on the Euler-Bernoulli beam theory (see [2,12,20,23]), while the last term results from the geometric nonlinearity included in Eq. (1). We restrict our attention to the first mode of vibration, hence the deflection is approximated by where w a (t) is the modal amplitude, and ϕ(x) is the mode shape function assumed to have the form Such a mode shape allows to incorporate (indirectly) into the model the effect of the end moments M 2 . The parameter α 0 will be further specified as a function of other characteristics of the system. For the triangular frame member α 0 > 0, while the special case α 0 = 0 corresponds to a simply supported beam (see Fig. 3).
Using the classical Galerkin method, we obtain the approximate equation of motion for the modal amplitude: where q = w a /d 1 and N 2 is the relative (dimensionless) axial force. The coefficients c 1 and c 3 are dependent on the beam parameters, and can be determined by considering the static case (q = 0) and the relationships given in [14]: where Note that N E 2cr is the Euler critical axial load for a beam/column, while N 2cr is the critical force that corresponds to the static critical load F cr calculated for the triangular frame. The quantity N 2 in Eq. (10) is normalized just by N 2cr , that is, In the examples below, we suppose that N 2 (t) is a periodic function: where T * = 2π/ω. It is convenient to introduce the dimensionless quantities: where ω 1 = √ c 1 is the first natural frequency of a load-free beam. Finally, Eq. (10) can be rewritten as Now q = q(τ ), an overdot denotes differentiation with respect to τ , and N 2 (τ ) is periodic with a period T = ω 1 T * = 2π/Ω. Thus, the fully nondimensionalized ordinary differential equation (15) takes a form of the well-known Hill's equation enriched with a nonlinear term of the Duffing type. This discrete system can be treated as a minimal model for vibrations of the triangular frame, focused on its neuralgic member, i.e. the cross-beam.

Remarks on analytical treatment of nonlinear time-periodic ODEs
The single approximate governing equation (15) seems to have a quite simple form. However, it combines two fundamental features that considerably complicate solving the problem analytically. Firstly, the model is an ODE with a time-periodic coefficient, i.e. a parametrically excited system. Secondly, it includes a strong nonlinearity. For the last decades, especially since 1950s, extensive efforts have been undertaken in the area of linear parametric systems, where the Floquet theory plays the central role. In particular, due to its usefullness in many physical problems, the Mathieu equationq has been studied thoroughly and elucidated. Various methods have been employed to investigate stability of steady-state solutions. The stability boundaries can be approximated by using, among others, i) infinite Hill's determinant, the harmonic balance method, or ii) perturbation methods (e.g. the Lindstedt-Poincaré technique, the method of multiple scales) [2][3][4][5]8,17,18]. The second approach requires the assumption that the time-varying coefficient is relatively small (ε 1). In the context of beams under axial excitation, the Mathieu-Hill type equations are analysed, for instance, in [9,10,23].
Today there exist a number of approximate analytical solution procedures applicable to nonlinear dynamical systems (see [1,6,18]). Nevertheless, the vast majority of them take their roots from the classical perturbation approach (e.g. the averaging methods, the Krylov-Bogolubov-Mitropolskij method, the method of multiple  scales, the homotopy perturbation method, etc.). Thus, the techniques work well in case of weakly nonlinear oscillatorsq when the so-called small parameter ε 1 can be associated with the nonlinear part Q of the governing equation. In recent years, some analytical methods for strongly nonlinear systems have been developed, too (e.g. the modified/adopted Lindstedt-Poincaré method or homotopy method). Usually, they are suited to certain type of nonlinearity only and are not generally applicable [6,16].
As will be seen further, from practical viewpoint Eq. (15) cannot be considered to contain the small perturbation parameter. In other words, the influence of both the cubic and time-periodic terms is rather strong. Due to limitations of the most of analytical methods and for convenience reasons, character of our studies is almost purely numerical.

Vibrations under harmonic load
We focus attention on a symmetrical triangular frame with the parameter values given in Table 1. This particular set of values is used in the railway industry. The area of the arms' cross section, A 1 , is constant, while the dimension b (and c = A 1 /b) will be altered.
In engineering practice, the width of the arms is one of the most important parameters. In the static case, for example, flat buckling of the triangular frame can occur for small values of b (at the given data: b ≤ 18 mm); otherwise, the system undergoes lateral buckling [14]. As presented in Fig. 4, the static critical load of the beam, N 2cr , increases with increasing b. At the same time, such a strengthening of the arms in respect of in-plane bending decreases the coefficient of the cubic term, γ . The question is how it affects the dynamics of the beam.
First, we assume that the structure is subjected to harmonic load and the relative axial force on the crossbeam is N 2 (t) = a sin(ωt) with ω = 2π f . Consequently, the initial value problem under consideration is The quantity ε 0 can be regarded as the relative amplitude of the beam imperfection. Initially, the results obtained for b = 20 mm are presented (γ ≈ 1.69); suppose that the imperfection amplitude ε 0 = 10 −3 . Figure 5a shows the parameter plane (Ω, a) with regions of stable and unstable motion of the system. More precisely, the areas filled with grey relate to positive values of the maximal Lyapunov exponent (MLE), λ 1 . The map has been computed with steps f = 1 Hz and a = 0.1. The Lyapunov exponents have been evaluated using the Gram-Schmidt orthogonalization procedure and renormalization [13,19]. Taking into account relatively slow convergence of the numerical computations in the analysed undamped case, vibrations are classified as chaotic if λ 1 ≥ ε 1 , where ε 1 is a small positive value (here ε 1 = 0.02).
For comparison purposes, stability boundaries for the abbreviated system (γ = 0) are shown. The transition curves correspond to 1T -periodic (solid line) and 2T -periodic solutions (dashed line) [4,8,17]. Here, the dependencies a(Ω) have been determined by means of the Rayleigh method that, in fact, is a variant of the harmonic balance method. The steady-state solutions are assumed to have the forms [8]: and As can be seen, the numerical results based on the nonlinear model fit well into the instability tongues for the linear system. However, these two cases are not fully comparable. Unstable solutions of the Hill equation (or the Mathieu equation) increase unboundedly with time, while the nonlinear term γ q 3 causes amplitude limitation [17]. Thus, deflections of the beam may become very large but remain bounded. This type of behaviour can be detected by analysis of the displacement and velocity amplitudes. Let X = [q,q] T be the state vector of the system. Character of the system evolution can be assessed, for example, by the measure where · denotes the Euclidean norm. The maximal values of X rel (τ ) found in the range 0 ≤ τ ≤ 5 × 10 4 are plotted in Fig. 5b. The contour line corresponding to X rel = 100 coincides largely with the curves a(Ω), especially for lower a. As the load amplitude increases, X rel grows significantly but does not exceed 4 × 10 4 .  As results from Fig. 5a, there are scattered subregions of regular motion for a > 1, above the curves corresponding to the linear system. In fact, quasi-periodic and almost periodic vibrations can be easily found there.
By analogy to Fig. 5, a stability diagram for b = 50 mm (γ ≈ 1.14) is shown in Fig. 6. As in the previous case, the critical value of λ 1 has been set to 0.02. Again the lower contour of the unstable region is quite well approximated by the curves determined analytically for the linear system. Nevertheless, it seems that an increase in b gradually enlarges the scattered regions indicating regular dynamics and shifts them towards lower values of both the load amplitude and frequency.
More detailed analysis of numerical solutions is presented in the next subsection.

Vibrations under pulsating load
The example of harmonic load is convenient for confrontation of the results with the classical linear case. Nevertheless, the minimal model based on the relationships given in Ref. [14] is dedicated to a cross-beam being in compression rather than in tension. Therefore, let us assume now that the beam is subjected to a periodic compressive axial load, and the following initial value problem is considered: so that the period of the load is still T = 2π/Ω. In Fig. 7a, the stable and unstable regions for the beam are presented for b = 20 mm (λ 1 = 0.02, ε 0 = 10 −3 ). The transition curves resulting from the simple approximation of type (20) do not occur in pairs and therefore do not form the typical tongues. However, the stability boundary with a tip at Ω = 2 is sufficiently accurately predicted by the linear case. Unlike in the previous numerical experiment, there are no scattered subregions of stable solutions for high values of a.
To give a deeper insight into the approximate values of the maximal Lyapunov exponent, a more detailed filled contour plot is given in Fig. 7b. Each cell is a rectangle of dimensions Ω × a. As can be seen, there are only several single points representing λ 1 < 0. The stable region relates mostly to the range 0 ≤ λ 1 < 0.01. The unstable region, in turn, is mainly occupied by the cells connected with evidently chaotic dynamics (λ 1 ≥ 1). Thus, the maps should be treated as a rather rough evaluation tool for determination of the vibration type. Naturally, this approximate character is a matter of a balance between the three factors: resolution of the map ( Ω, a), length of the analysed time range of the system evolution, and computation time. Selected steady-state responses of the system for Ω = 1.6 (see Fig. 8) are illustrated in Fig. 9 in the form of phase portraits. Similar to the linear case, the stable region is dominated by aperiodic vibrations. Below and above the instability tongue (a = 0.3 and a = 1) quasi-periodic motion occurs (Fig. 9a, c). In the first case, the plane phase (q,q) is occupied by the trajectory projection so densely that the diagram constitutes an annulus. More disordered vibrations with much higher amplitudes of both the displacement and velocity can be observed within the unstable regions. The solution corresponding to a = 2 is depicted in Fig. 9b, while apparently chaotic behaviour for a = 3 is shown in Fig. 9d. Note that q > 1 means w a > d 1 , i.e. large deflections of the beam, which definitely goes beyond the classical (linear) beam theory.
The vibrations can be more precisely classified by means of the Poincaré maps and frequency spectra. The results of the application of the fast Fourier transform (FFT) are presented in Fig. 10. Here, Ω k is a discrete dimensionless frequency (k = 0, 1, 2, . . .), i.e. related to ω 1 , whereas |q k | denotes the spectral magnitude. The spectra of a discrete nature with two dominating frequencies correspond to the quasi-periodic solutions (Fig. 10a, c). Presumably, these frequencies are incommensurable. The chaotic solution typically has a continuous Fourier spectrum (Fig. 10d). A similar effect occurs for the other unstable case; however, the spectral lines are concentrated over a narrower frequency range (Fig. 10b).
Instability of the latter type may require more attention. As mentioned before, the beam is excited to high amplitudes which nevertheless do not increase exponentially like in the linear system. Differences between such behaviour and the stable one are shown in Fig. 11. In the standard quasi-periodic case, the peak amplitude varies approximately sinusoidally. The unstable response undergoes another kind of modulation: the amplitude changes sharply from very low to very high values (small vs. large deflections). This beatinglike phenomenon was analysed by Kidachi and Onogi [11] who studied the Mathieu equation modified by cubic nonlinearity. The researchers used the term 'pseudo-unstable state'. More precisely, such a solution is unstable in the sense that it exhibits exponential growth and decay by turns. For the same reason, the response does not escape to infinity, and the system cyclically returns to the neighbourhood of the initial state.
Comparison of the discussed results and the ones for b = 50 mm indicated very slight differences in the stability diagrams and the maps of the MLE values. Also, a more detailed analysis of the MLE computed at selected levels of the amplitude a in the frequency range 0 < Ω ≤ 2.5, for 10 ≤ b ≤ 90 mm ( b = 1 mm) has been conducted. The effect of b on λ 1 seems to be negligible; all in all, the stable and unstable ranges (at certain a) remain the same. Moreover, the results turned out to be almost insensitive to changes in the initial value, ε 0 . However, it should be noted that only relatively small imperfection amplitudes were examined (e.g. ε 0 = 2 × 10 −3 , 5 × 10 −3 , 10 −2 ).

Conclusions
The problem of dynamics of the symmetrical triangular frame has been reduced to the problem of transverse parametric vibration of the cross-beam. The minimal mathematical model of the system has been formulated by means of the Hamilton principle and the Galerkin method. The time-periodic and nonlinear character of the system has been briefly discussed in terms of approximate analytical methods. The regions of stable and unstable solutions have been determined in the 'load frequency-amplitude' plane. For harmonic axial load, the lower boundaries of the parametric instability zones are well approximated by the transition curves specified for the abbreviated system. Moreover, some scattered subregions of regular motion occur and become larger as the width of the frame arms increases. This effect is not present in the case of pulsating load. The stable regions are dominated by quasi-periodic dynamics of the nonlinear system. Within the unstable regions, in turn, the vibration amplitude of the beam increases significantly but remains bounded. The behaviour is not necessarily strictly chaotic, the pseudo-unstable states can be observed too.
This kind of analysis, focused on both the vibration type and amplitude growth, is crucial in the formulation of safety conditions for the structure, i.e. certain criteria for the so-called technical stability. However, from a practical viewpoint, the dynamic stability issue is associated mainly with behaviour of a system under suddenly applied loads (e.g. an impulse or step load). Since the problem is of great importance in the particular case of brake triangles, it merits special attention in further studies of the mechanical system.