Bending–torsional-coupled vibrations and buckling characteristics of single and double composite Timoshenko beams

Abstract

The stability and free vibration analyses of single and double composite Timoshenko beams have been investigated. The closed-section beams are subjected to constant axially compressive or tensile forces. The double beams are assumed to be connected by a layer of elastic translational and rotational springs. The coupled governing partial differential equations of motion are discretized, and the resulted eigenvalue problem is solved numerically by applying the Chebyshev spectral collocation method. The effects of the elastic layer parameters, the axial forces, the slenderness ratio, the bending–torsional coupling, and the boundary conditions on the critical buckling loads, mode shapes, and natural transverse frequencies have been studied. A parametric study was performed, and the obtained results revealed different features, which hopefully can be useful for single- and double-beam-like engineering structures.

Introduction

Beam-type structures are basic components of several applications in the mechanical and aerospace engineering disciplines and are commonly used in constructions and structural systems. Understanding the different characteristics of the structures is crucial to ensure a good performance and to avoid failure. For example, determining the natural frequencies of a structure is important to avoid resonance. Moreover, estimating the critical buckling loads is essential to prevent instability and to help the designer in choosing an appropriate length and cross-sectional area of the beam to be used in such application.

Due to their high strength-to-weight ratio, high stiffness, and vibration absorption; double-beam systems are widely used in engineering applications such as sandwich or composite beams, multiple-walled carbon nanotubes, aircraft wing spars, tall building, double-beam cranes, railway tracks resting on elastic foundations, pipelines, bridge spans, and trusses. ¹ In some applications, as the main beam is excited, passive vibration control could be achieved by suppressing the undesired vibration by the aid of the auxiliary (secondary) beam.^2,3 Besides, the relative motion of the two beams could play an important role in dissipating the mechanical energy of the system.

Many researchers addressed the dynamic analysis of double-beam systems. For example, H Fei et al. ⁴ investigated the dynamic characteristics of a inclined and tensioned double-beam system. The two different beams are connected by elastic springs, and the governing equations of motion are derived by considering the effects of the flexural rigidity, sag, inclined angle, and the boundary conditions. The system was simplified into 4-degrees-of-freedom (DOF) system, and a semi-theoretical semi-numerical approach was proposed to solve the dynamical properties and the frequency equation of the structure. Furthermore, Murmu and Adhikari ⁵ presented the vibration analyses of coupled nanobeam system under initial compressive pre-stressed condition. Eringen’s nonlocal elasticity model was used to obtain the expressions for bending-vibration of the double-nanobeam-system. The nano-scale and coupling spring effects were examined, and it was shown that they have a significant influence on the natural frequencies of the nonlocal double-nanobeam-system. The results revealed that the nonlocal effects and the increase in the pre-stressed load tend to decrease the frequencies. In another study, Szmidt ⁶ experimentally and analytically carried out the dynamics of twin cantilever beams connected by a viscoelastic member, and the shear deformation was taken into consideration. A continuous system was introduced to identify the model by finding the stiffness and the damping properties of the free vibrations of a single cantilever beam excited in the first mode. The effect of the viscoelastic member on the vibration characteristics of the system was examined. Li et al. ⁷ discussed the free and forced vibration analysis of a double-beam system connected by a viscoelastic layer. A semi-analytical method was applied to find the natural frequencies and mode shapes. The forced vibration responses were obtained by utilizing a modal-expansion iterated method by deriving an orthogonality condition for the double-beam system to decouple the governing differential equations of motion. Moreover, KV Nguyen ⁸ studied the effect of a concentrated mass position on the free vibrations of a cracked double-beam system connected by an elastic layer. Numerical simulations were conducted, and it was shown that the frequency of the cracked double beam varies irregularly when the concentrated mass is located at the crack position, and the crack position can be detected by the location of peaks of the wavelet transform of the relationship between the natural frequency and the location of concentrated mass. Kelly and Srinivas ⁹ applied a normal mode solution to determine the natural frequencies and mode shapes for a set of Euler–Bernoulli beams subjected to axial loads and connected by continuous elastic layer. Mirzabeigy et al. ³ utilized the differential transform method to obtain the natural frequencies and the mode shapes of two parallel beams connected by a layer with a linear and parabolic variation. It was assumed that the ends of the beams were supported by translational and rotational springs. The influence of the stiffness layer, flexural rigidity of the beams, and the boundary conditions on the vibration characteristics was addressed. Oniszczuk ¹⁰ applied the classical Bernoulli–Fourier method to determine the natural frequencies of two parallel simply supported beams continuously joined by a Winkler elastic layer. Two partial differential equations (PDEs) were used to describe the motion of the system, and the initial-value problem was considered to find the natural frequencies and the mode shapes. In another study, Di Lorenzo et al. ¹¹ carried out the study of flexural behavior of layered elastically connected Euler–Bernoulli beams. It was assumed that the beams were supported with translational and rotational springs (considered as supports and joints). The theory of generalized function was used to derive the exact beam modes for the systems with discontinuities associated with the supports and joints. By applying the pertinent orthogonality conditions, the time response was found using the linear superposition and the modal analyses.

Moreover, the effect of the compressive force on the dynamics of beams undergoing bending–torsional-coupled vibrations was investigated by Prokić and Lukić. ¹² A closed-form solution of the equations of motion simply supported thin-walled beams was derived using Vlasov theory. The impact of axial force, bending–torsion coupling, rotary inertia, and warping stiffness on the behavior of thin-walled beams was demonstrated. In addition, Li et al. ¹³ applied the dynamic transfer matrix method to analyze the free vibration behavior of thin-walled Timoshenko beams subjected to axial loads. A continuous model was assumed, and the equations of motion for coupled bending and torsional vibration were obtained. A parametric study was carried out to determine the impacts of the axial compression force, warping stiffness, shear deformation, and rotary inertia on the frequencies and mode shapes of the Timoshenko beams.

In another study, J Li et al. ¹⁴ utilized the dynamic transfer matrix method to examine the free bending–torsional vibrations of axially loaded Euler–Bernoulli beams with thin-walled open monosymmetrical cross sections. The coupled vibration characteristics of axially loaded slender composite beams were investigated by J Li and X Jin. ¹⁵ The natural frequencies and mode shapes are obtained by solving the governing differential equations of motion. The normal mode method was used to obtain the displacement response of slender composite beams with solid and closed thin-walled cross sections.

Other articles addressed the mechanics of bending–torsional behavior of beams. For example, Vo and Lee ¹⁶ utilized the stationary value of total potential energy to examine the flexural behavior of composite box beams. The finite element method was applied to predict the deflection of thin-walled composite box beams acted upon by vertical and torsional loads. The effects of boundary conditions, laminate orientation, and span-to-height ratio of the composite beam on the vertical displacement and angle of twist were illustrated.

Also, Latalski et al. ¹⁷ discussed the free and forced vibrations of a rotating rigid hub and a flexible composite thin-walled beam. Hamilton’s principle is used to derive the PDEs of motion of the beam that is featuring bending-twist coupling. The influence of the hub mass, rotary inertia, and fiber lamination angle on the natural frequencies was considered. Sousa and Da Silva ¹⁸ carried out the behavior of multilayered composite beams with horizontal interlayer slip, where the displacements and slips were obtained after applying the boundary conditions. Kahya ¹⁹ presented the buckling response of laminated composite beams. The first-order shear deformation theory was utilized, and the finite element method was applied to obtain the critical buckling loads.

Based on the above review, it is noted that the free vibrations and the buckling characteristics of double composite Timoshenko beam systems were not investigated. In addition, in most cases, it was assumed that the elastic layer, that connects the beams, is a set of translational springs only. In this theoretical study, the free vibration and buckling behavior of double composite Timoshenko beam systems, connected with an elastic layer that consists of translational and rotational springs, are proposed. The double-beam system is featuring bending–torsion coupling, which is of particular interest from an aeroelastic point of view. The spectral Chebyshev collocation method is used to discretize the governing PDEs of motion and to convert them into a system of algebraic equations. Then, the resulted eigenvalue problem is solved to determine the natural frequencies, mode shapes, and the critical buckling loads of the double-beam systems under consideration.

Mathematical model of a single and double composite Timoshenko beams

In order to obtain the governing equations of motion of a single straight beam subjected to axial compressive force, and featuring coupled bending–torsional vibrations, Hamilton’s principle will be utilized. The total potential energy, U and the kinetic energy, T of the axially loaded composite Timoshenko beam are given as

$$\begin{array}{c}\multicolumn{1}{c}{U=\frac{1}{2}\underset{0}{\overset{L}{\int }}\{E{I}_y{(\theta \prime )}^2-P[{(w\prime )}^2+(I_s/\mu ){(\psi \prime )}^2]}\end{array}+2K\theta \prime \psi \prime +\kappa \mathrm{AG}{(w\prime -\theta )}^2+\mathrm{GJ}{(\psi \prime )}^2\}\mathrm{dx}$$

(1a)

$$T=\frac{1}{2}\underset{0}{\overset{L}{\int }}\left\{\mu {\left(\stackrel{.}{w}\right)}^2+I_s{\left(\stackrel{.}{\psi }\right)}^2+\rho I_y{\left(\stackrel{.}{\theta }\right)}^2\right\}dx$$

(1b)

where ρ is the density of the material, E is Young’s modulus, $I_y$ is the second moment of area of the beams cross section about the y-axis, κ is the shear correction factor, A is the cross-sectional area, G is the shear modulus, K is the flexure–torsion coupling rigidity, $\mu =\rho A$ is the mass per unit length, P is the constant axial force passes through the centroid of the cross section, $I_s$ is the polar mass moment of inertia per unit length about the x-axis, J is the polar moment of area, θ is the bending rotation about the z-axis, w is the flexural displacement in the z-direction, and ψ is the torsional rotation about the x-axis. Thus, EI_y, GJ, K, and κGA are the flexural rigidity, torsional rigidity, flexure–torsion coupling rigidity, and shear rigidity of the composite beam, respectively. ²⁰ The dot and the prime denote the differentiation with respect to time t and the spatial axial coordinate x, respectively. In the analysis, the constant axial force P is considered to be positive when it is compressive as shown in Figure 1 and considered to be negative when it is tension. ²⁰

Figure 1.

The double composite Timoshenko beam system under consideration.

Hamilton’s principle is applied to obtain the governing equations of motion and the boundary conditions

$$\underset{t_1}{\overset{t_2}{\int }}(\delta T-\delta U+\delta W)\mathrm{dt}=0$$

(2)

where W is the virtual work done by the nonconservative forces, which is set to zero in this study. Substituting equations (1a) and (1b) into equation (2), and setting $\delta w=\delta \psi =\delta \theta =0$ at $t_1$ and $t_2$ , the governing PDEs of motion of a straight composite Timoshenko beam (with length L, height h, and breadth b) subjected to an axial force P and featuring bending–torsion coupling are given as ²⁰

$$-\rho I_y\stackrel{··}{\theta }+E{I}_y\theta ″+\kappa \mathrm{AG}(w\prime -\theta )+K\psi ″=0$$

(3a)

$$-\mu \stackrel{··}{w}-Pw″+\kappa \mathrm{AG}(w″-\theta \prime )=0$$

(3b)

$$-I_s\stackrel{··}{\psi }-\frac{P{I}_s}{\mu }\psi ″+K\theta ″+\mathrm{GJ}\psi ″=0$$

(3c)

Due to the symmetry in the xz and xy planes, and as the x-axis is the locus of the shear centers, only the ply orientation causes the bending–torsional vibrational coupling. The equations were derived previously in Kaya and Ozgumus ²⁰ and repeated in this study for convenience. Assuming harmonic solutions in time as

$$w(x,t)=\overline{w}\left(x\right)e^{i\omega t}$$

(4a)

$$\theta (x,t)=\overline{\theta }\left(x\right)e^{i\omega t}$$

(4b)

$$\psi (x,t)=\overline{\psi }\left(x\right)e^{i\omega t}$$

(4c)

where $i=\sqrt{-1}$ , and ω is the circular natural frequency. The following dimensionless parameters are introduced as ²⁰

$$W=\frac{w}{L},X=\frac{x}{L},r^2=\frac{I}{A{L}^2}$$

(5)

Substituting equations (4a)–(4c) and (5) into equation (3a)–(3c) yields

$$\frac{E{I}_y}{L^2}\theta ″+\kappa \mathrm{AG}(W\prime -\theta )+\frac{K}{L^2}\psi ″=-\rho I_y{\omega }^2\theta$$

(6a)

$$-\frac{P}{L}W″+\frac{\kappa \mathrm{AG}}{L}(W″-\theta \prime )=-\rho \mathrm{AL}{\omega }^{2}W$$

(6b)

$$-\frac{P{I}_s}{\mu L^2}\psi ″+\frac{K}{L^2}\theta ″+\frac{\mathrm{GJ}}{L^2}\psi ″=-I_s{\omega }^2\psi$$

(6c)

In this study, two types of boundary conditions, free (F) and clamped (C) are considered. For a clamped edge, the displacements W, θ, and ψ equal zero; whereas for a free edge, the bending moment, shear force, and the torque equal to zero, and they are given as: ²⁰

• Bending moment

$$E{I}_y\theta \prime +K\psi \prime =0$$

(7a)

• Shear force

$$-\mathrm{Pw}\prime +\kappa \mathrm{AG}(w\prime -\theta )=0$$

(7b)

• Torque

$$-(P{I}_s/\mu )\psi \prime +K\theta \prime +\mathrm{GJ}\psi \prime =0$$

(7c)

For the double-beam system, equations (4a) and (4b) can be applied to each of the beams that are connected by a continuous elastic layer (that consists of linear translational and rotational springs) shown in Figure 1. Assuming that the beams are identical (have the same material properties), the governing equations of motion are given as

$$\begin{array}{c}\multicolumn{1}{c}{\frac{E{I}_y}{L^2}{\theta }_1^{″}+\kappa \mathrm{AG}({W\prime }_1-{\theta }_1)+\frac{K}{L^2}{\psi }_1^{″}+K_{\theta }({\theta }_1-{\theta }_2)}\\ \multicolumn{1}{c}{=-\rho I_y{\omega }^2{\theta }_1}\end{array}$$

(8a)

$$\begin{array}{c}\multicolumn{1}{c}{-\frac{P}{L}{W}_1^{″}+\frac{\kappa \mathrm{AG}}{L}(W″-{\theta \prime }_1)+K_{W}L(W_1-W_2)}\\ \multicolumn{1}{c}{=-\rho \mathrm{AL}{\omega }^2{W}_1}\end{array}$$

(8b)

$$-\frac{P{I}_s}{\mu L^2}{\psi }_1^{″}+\frac{K}{L^2}{\theta }_1^{″}+\frac{\mathrm{GJ}}{L^2}{\psi }_1^{″}+K_{\psi }({\psi }_1-{\psi }_2)=-I_s{\omega }^2{\psi }_1$$

(8c)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{E{I}_y}{L^2}{\theta }_2^{″}+\kappa \mathrm{AG}({W\prime }_2-{\theta }_2)+\frac{K}{L^2}{\psi }_2^{″}+K_{\theta }({\theta }_2-{\theta }_1)}\\ \multicolumn{1}{c}{=-\rho I_y{\omega }^2{\theta }_2}\end{array}$$

(8d)

$$\begin{array}{c}\multicolumn{1}{c}{-\frac{P}{L}{W}_2^{″}+\frac{\kappa \mathrm{AG}}{L}({W″}_2-{\theta \prime }_2)+K_{W}L(W_2-W_1)}\\ \multicolumn{1}{c}{=-\rho \mathrm{AL}{\omega }^2{W}_2}\end{array}$$

(8e)

$$-\frac{P{I}_s}{\mu L^2}{\psi }_2^{″}+\frac{K}{L^2}{\theta }_2^{″}+\frac{\mathrm{GJ}}{L^2}{\psi }_2^{″}+K_{\psi }({\psi }_2-{\psi }_1)=-I_s{\omega }^2{\psi }_2$$

(8f)

where K_θ and K_ψ are the stiffnesses of the rotational springs that resist the rotation around the z and x axes, respectively, and K_W is the translational stiffness of the spring that resists the translation in the z-axis. It is assumed that the stiffnesses of the translational/rotational springs are independent of the material properties of the beams.

Solution procedure: Chebyshev collocation method

When performing analysis on single and double composite Timoshenko beam undergoing coupled bending–torsion coupling, it becomes challenging to obtain a closed-form solution to the vibration and buckling problems for these systems. Hence, the Chebyshev spectral collocation method is utilized for numerical solutions. Gauss–Chebyshev–Lobatto or Chebyshev collocation points are defined as the projections of equally spaced segments on the upper half of the unit circle on the interval [–1, 1]. These points are numbered from right to left, and are defined as ²¹

$$x_j=\cos (j\pi /N),j=0,1,\dots ,N$$

(9)

The Chebyshev differentiation matrix, [D]_N , is of size $(N+1)\times (N+1)$ and is formulated by interpolating a Lagrange polynomial of degree N at each Chebyshev point, differentiating the polynomial, and then evaluating the polynomial at each Chebyshev point. The elements of the [D]_N matrix are given as ²¹

$$\begin{array}{c}{\left(D_N\right)}_{00}=\frac{2N^2+1}{6},{\left(D_N\right)}_{\mathrm{NN}}=-\frac{2N^2+1}{6},\hfill \\ \multicolumn{1}{c}{{\left(D_N\right)}_{\mathrm{jj}}=\frac{-x_j}{2(1-x_j^2)},j=1,\dots ,N-1}\\ {\left(D_N\right)}_{\mathrm{ij}}=\frac{c_i}{c_j}\frac{{(-1)}^{i+j}}{(x_i-x_j)},i\ne j,i,j=0,\dots ,N.\hfill \\ \multicolumn{1}{c}{c_i=\{\begin{array}{cc}2,\hfill & i=0\mathrm{or}N\hfill \\ \multicolumn{1}{c}{1,}& \mathrm{otherwise}\hfill \end{array}}\end{array}$$

(10)

In this study, the X axis is normalized to be in the range of [0, 1]; thus, the Chebyshev differentiation matrices will have different entries than those obtained by Trefethen21 as they depend on the distribution of the points. Due to its accuracy, fast convergence, and simplicity in implementation, the Chebyshev spectral collocation method has effectively been utilized to carry out the free vibration and buckling analyses of one, two, and three-dimensional continuous systems (at the macro- and micro-levels).^22–24 When discretizing ordinary differential equations (ODEs) and PDEs by the Chebyshev collocation method, the nth derivative of a function is represented by Dn = (D_N) ⁿ . According to the Timoshenko beam model under consideration, each Chebyshev point possesses 3 DOFs—θ, Ψ, and W. Hence, a displacement vector is defined as

$$\begin{array}{c}\multicolumn{1}{c}{[U]=[{\theta }_1{\psi }_1{W}_1{\theta }_{N+1}{\psi }_{N+1}{W}_{N+1}{\theta }_2\cdots {\theta }_N}\\ \multicolumn{1}{c}{{\psi }_2\cdots {\psi }_N{W}_2\cdots W_N{]}^T}\end{array}$$

(11)

where N + 1 is the number of points in the X direction (along the beam’s span); θ₁, Ψ₁, and W₁ are the displacements at X = 0, and θ_{N  +  1}, Ψ_{N  +  1}, and W_{N  +  1} are the displacements at X = 1. The size of the vector {U} is $3M\times 1$ , where M = N + 1. The left-hand sides of equations (6a)–(6c) are discretized by the Chebyshev collocation method as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}1=(\frac{E{I}_y}{L^2}*\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\kappa \mathrm{AG}*\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes D1)}\\ \multicolumn{1}{c}{-(\kappa \mathrm{AG}*\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\end{array}$$

(12a)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}2=(\frac{-P}{L}*\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{\kappa \mathrm{AG}}{L}*\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{-(\frac{\kappa \mathrm{AG}}{L}*\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D1)}\end{array}$$

(12b)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}3=(\frac{-P{I}_s}{\mu L^2}*\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{\mathrm{GJ}}{L^2}*\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\end{array}$$

(12c)

In a similar manner, the right-hand sides of equations (6a)–(6c) are discretized as

$$\mathrm{RH}1=-\rho I_y*\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I$$

(13a)

$$\mathrm{RH}2=-\rho \mathrm{AL}*\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes I$$

(13b)

$$\mathrm{RH}1=-I_s*\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes I$$

(13c)

The whole system of equations is put in matrix vector form as

$$\left[\begin{array}{c}\mathrm{GE}1\hfill \\ \multicolumn{1}{c}{\mathrm{GE}2}\\ \multicolumn{1}{c}{\mathrm{GE}3}\end{array}\right]\left\{\begin{array}{c}\theta \hfill \\ \multicolumn{1}{c}{\psi }\\ \multicolumn{1}{c}{W}\end{array}\right\}={\mathrm{\Omega }}^2\left[\begin{array}{c}\mathrm{RH}1\hfill \\ \multicolumn{1}{c}{\mathrm{RH}2}\\ \multicolumn{1}{c}{\mathrm{RH}3}\end{array}\right]\left\{\begin{array}{c}\theta \hfill \\ \multicolumn{1}{c}{\psi }\\ \multicolumn{1}{c}{W}\end{array}\right\}$$

(14)

The size of each term (GE1, GE2, GE3, RH1, RH2, and RH3) is $N\times 3N$ , and the size of the system given in equation (14) is $3N\times 3N$ . The boundary conditions (in dimensionless form) at the edge X = 0 are applied using the Chebyshev collocation method and the Kronecker product as follows:

• Clamped (C)

$$(\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes \left[\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right])[U^T]=0$$

(15a)

$$(\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes \left[\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right])[U^T]=0$$

(15b)

$$(\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes \left[\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right])[U^T]=0$$

(15c)

• Free (F)

$$\left(\right(\frac{E{I}_y}{L}*\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D1(1,:))+(\frac{K}{L}*[010]\otimes D1(1,:)))\left\{U\right\}=0$$

(16a)

$$\left(\left(\left(-P+\kappa AG\right)*\left[\begin{array}{ccc}0& 0& 1\end{array}\right]\otimes D1\left(1,:\right)\right)-\left(\kappa AG*\left[\begin{array}{ccc}1& 0& 0\end{array}\right]\otimes \left[\begin{array}{ccc}1& \begin{array}{cc}0& \cdots \end{array}& 0\end{array}\right]\right)\right)\left\{U\right\}=0$$

(16b)

$$\left(\left(\left(\frac{-P{I}_s}{\mu L}+\frac{GJ}{L}\right)*\left[\begin{array}{ccc}0& 1& 0\end{array}\right]\otimes D1\left(1,:\right)\right)+\left(\frac{K}{L}*\left[\begin{array}{ccc}0& 1& 0\end{array}\right]\otimes D1\left(1,:\right)\right)\right)\left\{U\right\}=0$$

(16c)

where $[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}]$ is a $1\times N+1$ vector, with the first entry equal to one and other entries equal to zero, since it corresponds to the boundary at X = 0, and D1(1,:) is the first row of the D1 matrix. When the boundary conditions are evaluated at X = 1, the vector $[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}]$ is replaced by the vector $[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}]$ whose size is $1\times N+1$ , and the last entry equals one and the others are zero, and the vector D1(1,:) is replaced by the vector D1(N + 1,:), which is the last row of the D1 matrix.

For the double-beam system, the displacement vector is defined as

$$\begin{array}{c}\multicolumn{1}{c}{\{U\}=\{{\theta }_{1,1}{\psi }_{1,1}{W}_{1,1}{\theta }_{2,1}{\psi }_{2,1}{W}_{2,1}{\theta }_{1,N+1}{\psi }_{1,N+1}}\\ \multicolumn{1}{c}{W_{1,N+1}{\theta }_{2,N+1}{\psi }_{2,N+1}{W}_{2,N+1}}\\ \multicolumn{1}{c}{{\theta }_{1,2}\cdots {\theta }_{1,N}{\psi }_{1,2}\cdots {\psi }_{1,N}{W}_{1,2}\cdots }\\ \multicolumn{1}{c}{W_{1,N}{\theta }_{2,2}\cdots {\theta }_{2,N}{\psi }_{2,2}\cdots {\psi }_{2,N}{W}_{2,2}}\\ \multicolumn{1}{c}{\cdots W_{2,N}{\}}^T}\end{array}$$

(17)

where the first subscript refers to the ith beam (beam 1 or beam 2), and the second subscript refers to the jth Chebyshev point along the beam’s span (j = 1, 2, …, N + 1). The governing equations (8a)–(8f) are discretized as

$$\begin{array}{c}GE1=\left(\frac{E{I}_y}{L^2}*\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\end{array}\right]\otimes D2\right)+\left(\kappa AG*\left[\begin{array}{cccccc}0& 0& 1& 0& 0& 0\end{array}\right]\otimes D1\right)\\ -\left(\kappa AG*\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\end{array}\right]\otimes I\right)+\left(\frac{K}{L^2}*\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\end{array}\right]\otimes D2\right)\\ +\left(K_{\theta }*\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\end{array}\right]\otimes I\right)-\left(K_{\theta }*\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\end{array}\right]\otimes I\right)\end{array}$$

(18a)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}2=(\frac{-P}{L}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{\kappa \mathrm{AG}}{L}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{-(\frac{\kappa \mathrm{AG}}{L}*\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D1)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{cccccc}0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(K_{W}L*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{-(K_{W}L*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes I)}\end{array}$$

(18b)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}3=(\frac{-P{I}_s}{\mu L^2}*\left[\begin{array}{cccccc}0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{\mathrm{GJ}}{L^2}*\left[\begin{array}{cccccc}0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(K_{\psi }*\left[\begin{array}{cccccc}0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{-(K_{\psi }*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes I)}\end{array}$$

(18c)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}4=(\frac{E{I}_y}{L^2}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\kappa \mathrm{AG}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes D1)}\\ \multicolumn{1}{c}{-(\kappa \mathrm{AG}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(K_{\theta }*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{-(K_{\theta }*\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\end{array}$$

(18d)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}5=(\frac{-P}{L}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{\kappa \mathrm{AG}}{L}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{-(\frac{\kappa \mathrm{AG}}{L}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D1)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(K_{W}L*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{-(K_{W}L*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\end{array}$$

(18e)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{GE}6=(\frac{-P{I}_s}{\mu L^2}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{K}{L^2}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(\frac{\mathrm{GJ}}{L^2}*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes D2)}\\ \multicolumn{1}{c}{+(K_{\psi }*\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\otimes I)}\\ \multicolumn{1}{c}{-(K_{\psi }*\left[\begin{array}{cccccc}0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\otimes I)}\end{array}$$

(18f)

The right-hand side of the governing equations and the boundary conditions are discretized following the same scheme. Once the boundary conditions are applied, the system (single or double beam) is written as

$$\left[\begin{array}{cc}\left[S_{\mathrm{BB}}\right]\hfill & \left[S_{\mathrm{BI}}\right]\hfill \\ \multicolumn{1}{c}{\left[S_{\mathrm{IB}}\right]}& \left[S_{\mathrm{II}}\right]\hfill \end{array}\right]\left\{\begin{array}{c}\left\{U_B\right\}\hfill \\ \multicolumn{1}{c}{\left\{U_I\right\}}\end{array}\right\}={\mathrm{\Omega }}^2\left[\begin{array}{cc}\left[0\right]\hfill & \left[0\right]\hfill \\ \multicolumn{1}{c}{\left[0\right]}& \left[Q\right]\hfill \end{array}\right]\left\{\begin{array}{c}\left\{U_B\right\}\hfill \\ \multicolumn{1}{c}{\left\{U_I\right\}}\end{array}\right\}$$

(19)

where the subscripts B and I indicate the grid points used for writing the collocation analog of the boundary conditions and the governing differential equation, respectively. For the single beam, the size of $S_{\mathrm{BB}}$ is $6\times (3M-6)$ , the size of S_BI is $(3M-6)\times 6$ , the size of S_IB is $(3M-6)\times (3M-6)$ , and the size of S_II is $(3M-6)\times (3M-6)$ . For the double-beam system the sizes of S_BB, S_BI, S_IB, and S_II are doubled. From the top partition of equation (19), we obtain

$$\begin{array}{c}\left[S_{\mathrm{BB}}\right]\left\{U_B\right\}+\left[S_{\mathrm{BI}}\right]\left\{U_I\right\}=\left[0\right]\hfill \\ \Rightarrow \left\{U_B\right\}=-{\left[S_{\mathrm{BB}}\right]}^{-1}\left[S_{\mathrm{BI}}\right]\left\{U_I\right\}\hfill \end{array}$$

(20)

Substituting equation (20) into equation (19) yields the final eigenvalue equation as

$$\begin{array}{c}\left[S_{\mathrm{IB}}\right]\left\{U_B\right\}+\left[S_{\mathrm{II}}\right]\left\{U_I\right\}={\mathrm{\Omega }}^2\left[Q\right]\left\{U_I\right\}\hfill \\ \Rightarrow (-\left[S_{\mathrm{IB}}\right]{\left[S_{\mathrm{BB}}\right]}^{-1}\left[S_{\mathrm{BI}}\right]+\left[S_{\mathrm{II}}\right])\left\{U_I\right\}={\mathrm{\Omega }}^2\left[Q\right]\left\{U_I\right\}\hfill \\ \Rightarrow {\left[Q\right]}^{-1}(-\left[S_{\mathrm{IB}}\right]{\left[S_{\mathrm{BB}}\right]}^{-1}\left[S_{\mathrm{BI}}\right]+\left[S_{\mathrm{II}}\right])\left\{U_I\right\}={\mathrm{\Omega }}^2\left\{U_I\right\}\hfill \end{array}$$

(21)

The natural frequencies are obtained by calculating the eigenvalues of the matrix

$$\left[A\right]={\left[Q\right]}^{-1}(-\left[S_{\mathrm{IB}}\right]{\left[S_{\mathrm{BB}}\right]}^{-1}\left[S_{\mathrm{BI}}\right]+\left[S_{\mathrm{II}}\right])$$

Results and discussion

A glass–epoxy composite beam with a rectangular cross section is used in the analysis, and it is assumed that the beam has unidirectional plies with fiber angles of 15°. The beam has the following properties (unless otherwise stated): ²⁰ width = 12.7 mm, thickness = 3.18 mm, EI = 0.2865 N m², GJ = 0.1891 N m², kAG = 6343.3 N, µ = 0.0544 kg/m, I_s = 0.777 × 10^–6 kg m, K = 0.1143 N m², L = 0.1905 m, r² = 0.00002322.

For the validation purpose, the natural frequencies of a single composite Timoshenko cantilever beam obtained from the proposed technique are compared with those obtained by Kaya and Ozgumus. ²⁰ The first six natural frequencies for clamped-free (cantilever) composite Timoshenko beams are computed using the presented numerical method where N is considered equal to 18 in the performed computations to attain convergence for the first six frequencies. As shown in Table 1, it is clear that the results are in good agreement with those found in Kaya and Ozgumus, ²⁰ which reflects the validity of the spectral collocation method. In addition, it is clear that the bending–torsion coupling tends to decrease the frequencies. It is worthwhile mentioning that in this study, the structure is considered as a closed box beam (beam with a closed cross section), and it is assumed that it has a negligible warping stiffness. Therefore, the bending–torsional vibrational coupling is only due to the ply orientation. However, for asymmetric composite beams with open sections (which is not the case in the present analysis), the coupling would be due to the warping of the cross section, the bending–torsion material coupling, and the non-coincidence of the elastic axis and the mass axis.^25,26

Table 1.

First six natural frequencies (Hz) of a single composite cantilever Timoshenko beam.

Natural frequencies
P = 7.5				P = –7.5
Coupled		Uncoupled		Coupled		Uncoupled
Currentstudy	Kaya andOzgumus 20	Currentstudy	Kaya and Ozgumus 20	Currentstudy	Kaya and Ozgumus 20	Currentstudy	Kaya andOzgumus 20
21.987	21.987	28.064	28.064	37.106	37.106	40.975	40.975
181.495	181.495	210.162	210.162	197.672	197.672	224.259	224.259
511.818	511.818	586.519	586.519	525.665	525.665	598.668	598.668
648.047	648.047	647.228	647.228	648.495	648.495	647.595	647.595
979.473	979.473	1113.952	1113.950	992.878	992.878	1125.711	1125.711
1558.134	1558.134	1766.662	1766.662	1571.338	1571.338	1778.429	1778.429

Figure 2 shows the effect of the square of the inverse of the slenderness ratio (r²) on the first three natural frequencies of a cantilever composite Timoshenko beam. It is observed that as the value of the inverse of the slenderness ratio increases, the natural frequencies decrease, where the decrease is more significant for the second and third modes. The figure reveals that the Euler–Bernoulli beam theory is sufficient for calculating the fundamental natural frequency for a thick closed box beam (as the inverse of the slenderness ratio has a negligible influence on the fundamental frequency), and the frequencies of a Timoshenko beam are less than those for an Euler–Bernoulli beam. This conclusion is true for closed box beams only with the properties given in Kaya and Ozgumus, ²⁰ whereas for thin-walled beams with closed cross sections, the discrepancies between Euler–Bernoulli and Timoshenko beam models may be significant even for very slender beams. ²⁷

Figure 2.

The effect of the inverse of the slenderness on the first three natural frequencies of a cantilever beam: (a) P = 7.5, K = 0; (b) P = 7.5, K = 0.1143; (c) P = –7.5, K = 0; and (d) P = –7.5, K = 0.1143.

The influence of the compressive axial force on the first three natural frequencies of a cantilever beam with r² = 0.01 is presented in Figure 3. It is observed that as the axial load increases, the natural frequencies decrease. In addition, one may notice that the bending–torsion coupling has a decreasing effect on the critical buckling loads of the beam. It is known that the critical buckling loads correspond to the zero natural frequencies. For example, the first three critical buckling loads of cantilever beam with the flexure–torsion coupling rigidity K = 0 are 19.5, 170.6, and 452.3 N compared to 14.8, 130, and 346.4 N for the cantilever beam with the flexure–torsion coupling rigidity K = 0.1143. The impact of the compressive axial force on the first three natural frequencies of a CC beam with r² = 0.01 is displayed in Figure 4, where CC denotes that the beam is clamped at both edges (X = 0 and X = 1). In this figure, the first three critical buckling loads with the flexure–torsion coupling rigidity K = 0 are 297.2, 574.3, and 1042 N compared to 226.8, 442, and 808.5 N for the CC beam with the flexure–torsion coupling rigidity K = 0.1143.

Figure 3.

Effect of axial force on the first three natural frequencies with the of a cantilever beam with r² = 0.01: (a) K = 0 and (b) K = 0.1143.

Figure 4.

Effect of axial force on the first three natural frequencies with the of a CC beam with r² = 0.01: (a) K = 0 and (b) K = 0.1143.

The first six non-normalized mode shapes of a cantilever composite Timoshenko beam subjected to an axial compressive load (P = 7.5) without and with the bending–torsion coupling are presented in Figures 5 and 6, respectively. From these figures, it can be concluded that the fourth mode shape is the fundamental torsion mode of the cantilever beams under consideration. Moreover, the first six mode shapes of the CC Timoshenko composite beam subjected to an axial compressive load (P = 7.5) with bending–torsion coupling, are presented in Figure 7.

Figure 5.

(a–f) First six non-normalized mode shapes of a composite cantilever Timoshenko beam (P = 7.5, r² = 0.00002322, K = 0).

Figure 6.

(a–f) First six non-normalized mode shapes of a composite cantilever Timoshenko beam (P = 7.5, r² = 0.00002322, K = 0.1143).

Figure 7.

(a–f) First six non-normalized mode shapes of a CC Timoshenko beam (P = 7.5, r² = 0.00002322, K = 0.1143).

Figure 8 is devoted to present the first three non-normalized mode shapes of CF-CF double-beam system exposed to an axial compressive load (P = 7.5), and the beams are connected by an elastic layer with K_W = 10 N/m² and K_Θ = K_Ψ = 10 N/rad. For these two identical beams (same cross-sectional area, same length, same material, and so on), it is clear that the first and third mode shapes (odd numbered modes) have the same magnitudes and opposite directions over the span of the beams. Also, the second mode shapes (even-numbered mode) of the beams have the same magnitudes and directions. The first three non-normalized mode shapes of a CF-CC double beams that are connected by an elastic layer with K_W = 100 N/m², and K_Θ = K_Ψ = 100 N/rad, are shown in Figure 9.

Figure 8.

(a–f) First three mode shapes (blue: θ, green: Ψ, red: W) of a CF-CF double beam, P = 7.5, K_W = 10 N/m², andK _$\Theta$  = K_Ψ = 10 N/rad; (a), (c), and (e) are for beam 1 and (b), (d), and (f) are for beam 2.

Figure 9.

(a–f) First three mode shapes (blue: θ, green: Ψ, red: W) of a CF-CC double beam, P = 0, K = 0.114, K_W = 100 N/m², and K _$\Theta$  = K_Ψ = 100 N/rad; (a), (c), and (e) are for beam 1 and (b), (d), and (f) are for beam 2.

In Figure 10, the effect of the axial compressive force on the first three natural frequencies of the CF-CF double-beam system with r² = 0.01, connected by an elastic layer with K_W = 10 N/m², and K_Θ = K_Ψ = 10 N/rad, is shown. It is noted that the natural frequencies decrease as the axial load increases. Moreover, the bending–torsion coupling has a decreasing effect on the critical buckling loads of the beam. For example, the first three critical buckling loads of a CF-CF double beam without coupling are 19.47, 170.6, and 274.3 N compared to 14.76, 31.75, and 129.9 N for the double beam with coupling. One may notice that the decrease in the critical buckling loads is more significant for higher modes. The effect of the axial compressive force on the first three natural frequencies of the CF-CF double beam at two different values of the inverse of the slenderness ratio (r²), and are connected by an elastic layer with K_W = 200 N/m², and K_Θ = K_Ψ = 200 N/rad, is presented in Figure 11. The figure reveals that at relatively high values of the elastic layer parameters, the value of r² has a negligible influence on the critical buckling loads.

Figure 10.

Effect of axial force on the first three natural frequencies of a CF-CF double composite Timoshenko beam with r² = 0.01, K_W = 100 N/m², and K _$\Theta$  = K_Ψ = 100 N/rad: (a) K = 0 and (b) K = 0.1143.

Figure 11.

Effect of axial force on the first three natural frequencies of a CF-CF double composite Timoshenko beam with K_W = 200 N/m², K _$\Theta$  = K_Ψ = 200 N/rad, and K = 0.1143: (a) r² = 0.00002322 and (b) r² = 0.01.

The effect of the axial force on the first three natural frequencies of the CC-CF double beam that are connected with an elastic layer with K_W = 10 N/m², and K_Θ = K_Ψ = 10 N/rad, is illustrated in Figure 12(a) and (b). In this case, the bending–torsion coupling has an increasing effect on the critical buckling loads. The first three critical buckling loads of the CF-CF double beam without coupling are 9.38, 160.7, and 288.2 N compared to 114.7, 210.4, and 333.5 N for the double beam with coupling. A point of interest is that the veering phenomenon occurs as the axial load increases. The veering refers to a region in which two eigenvalue loci approach each other and almost intersect as the system parameter is varied (the axial load in this case). Instead of intersecting, they veer away from each other. For instance, for the CC-CF double beam without coupling, at P = 28.14 N, there is veering between the second and the third natural frequencies. In addition, for CC-CF double beam with coupling, at P = 264.5 N, there is veering between the third and fourth natural frequencies. Furthermore, Figure 13(a) and (b) demonstrates the effect of the axial force on the first three natural frequencies of the CC-FF double beam with coupling at different parameters of the elastic layer. For the CC-FF beam connected by a layer with K_W = 10 N/m², and K_Θ = K_Ψ = 10 N/rad, the second and third critical buckling loads have the same values. The figure shows that as the stiffnesses of the springs of elastic layer increase, the critical buckling loads increase, and the veering phenomenon between the third and fourth natural frequencies take place at P = 100.9 N.

Figure 12.

Effect of axial force on the first three natural frequencies of a CC-CF double composite Timoshenko beam with K_W = 10 N/m², and K _$\Theta$  = K_Ψ = 10 N/rad, r² = 0.00002322: (a) K = 0 and (b) K = 0.1143.

Figure 13.

Effect of axial force on the first four natural frequencies of a CC-FF double composite Timoshenko beam with r² = 0.00002322, K = 0.1143: (a) K_W = K _$\Theta$  = K_Ψ = 10; (b) K_W = 100 N/m², and K_Ɵ = K_Ψ = 100 N/rad.

Tables 2 and 3 list the first six natural frequencies of the double-beam systems with different boundary conditions with and without coupling, respectively. The beams are connected by an elastic layer with K_W = 10 N/m², and K_Θ = K_Ψ = 10 N/rad, and are under the effect of an axial compressive load (P = 7.5). These tables reveal that the bending–torsion coupling has a decreasing or increasing effect on the natural frequencies of the presented double-beam systems. For example, the coupling has a decreasing effect on the CF-CF and CC-CC double-beam systems; whereas, it has an increasing effect on the FF-CF and CC-CF beams. The CF-CF double beam denotes that beam 1 is clamped at X = 0 and free at X = 1, and beam 2 is clamped at X = 0 and free at X = 1. For the beams with bending–torsion coupling, it is noted that the CF-CF beams have the lowest natural frequencies, while the CC-CC beams have the highest frequencies. Furthermore, the first, third, fifth, and sixth natural frequencies of the FF-CF double-beam system are higher than those of the CC-CF beams. However, the second and fourth natural frequencies of the FF-CF beam are lower than those for the CC-CF beams. From these tables, it is depicted that the CC-CF beams have the lowest natural frequencies and the CC-CC beams have the highest frequencies. Also, the first six frequencies of the FF-CF beams are lower than those for the CF-CF beams.

Table 2.

The first six frequencies of double Timoshenko beams (P = 7.5, r² = 0.00002322, K = 0.1143, K_W = 10 N/m², andK _$\Theta$  = K_Ψ = 10 N/rad).

Boundary conditions	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
CF-CF	4.424	21.987	158.583	181.495	487.841	511.818
FF-CF	130.929	180.023	475.313	487.825	515.591	957.256
CC-CF	121.963	185.841	220.492	492.300	509.924	954.666
CC-CC	177.244	188.463	494.293	508.861	952.029	966.089

Table 3.

The first six frequencies of double Timoshenko beams (P = 7.5, r² = 0.00002322, K = 0, K_W = 10 N/m², and K _$\Theta$  = K_Ψ = 10 N/rad).

Boundary conditions	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
CF-CF	28.064	190.264	210.162	571.711	586.519	647.228
FF-CF	19.743	187.093	209.977	478.149	572.745	589.834
CC-CF	11.291	178.979	198.617	214.275	571.260	582.186
CC-CC	209.188	215.841	570.182	578.436	1011.738	1082.731

In a similar manner, Tables 4 and 5 show the first six natural frequencies of double-beam systems with different boundary conditions with and without coupling, respectively. The beams are connected by an elastic layer with K_W = 10 N/m², and K_Θ = K_Ψ = 10 N/rad, and are subjected to an axial tensile load (P = –7.5). The coupling has a decreasing effect on the natural frequencies of the FF-CF and CC-CC double beams and an increasing effect on the natural frequencies of the CF-CF and CC-CF double beams.

Table 4.

The first six frequencies of double Timoshenko beams (P = –7.5, r² = 0.00002322, K = 0.1143, K_W = 10 N/m², and K _$\Theta$  = K_Ψ = 10 N/rad).

Boundary conditions	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
CF-CF	30.985	37.106	177.438	197.672	502.430	525.665
FF-CF	34.673	158.480	199.941	478.010	504.370	531.374
CC-CF	132.361	192.609	227.676	503.687	521.358	965.140
CC-CC	183.776	194.569	503.002	517.316	961.692	975.608

Table 5.

The first six frequencies of double Timoshenko beams (P = –7.5, r² = 0.00002322, K = 0, K_W = 10 N/m², and K _$\Theta$  = K_Ψ = 10 N/rad).

Boundary conditions	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
CF-CF	21.494	40.975	205.849	224.259	584.173	598.668
FF-CF	38.737	206.867	227.972	478.562	586.827	605.503
CC-CF	33.030	180.477	211.626	221.484	579.751	593.426
CC-CC	214.680	221.146	577.666	585.804	1012.676	1091.124

The frequencies for double-beam systems, which are not subjected to any axial force, are listed in Tables 6 and 7. It is seen that the coupling has a decreasing effect on the natural frequencies of the CF-CF, FF-CF, and CC-CC beams and an increasing effect on the CC-CF beams.

Table 6.

The first six frequencies of double Timoshenko beams (P = 0, r² = 0.00002322, K = 0.1143, K_W = 10 N/m², andK _$\Theta$  = K_Ψ = 10 N/rad).

Boundary conditions	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
CF-CF	22.623	30.747	168.295	189.779	495.187	518.791
FF-CF	24.138	145.497	190.158	476.982	496.034	523.386
CC-CF	127.364	189.281	223.991	498.143	515.561	959.953
CC-CC	180.543	191.543	498.668	513.108	956.873	970.861

Table 7.

The first six frequencies of double Timoshenko beams (P = 0, r² = 0.00002322, K = 0, K_W = 10 N/m², and K _$\Theta$  = K_Ψ = 10 N/rad).

Boundary conditions	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
CF-CF	35.283	198.214	217.341	577.976	592.626	647.411
FF-CF	30.725	197.441	219.003	478.356	579.916	597.642
CC-CF	24.988	179.729	205.520	217.644	575.606	587.752
CC-CC	211.954	218.511	573.937	582.133	1012.207	1086.936

Conclusion

The buckling behavior and the free vibration analysis of composite single and double Timoshenko beams, subjected to axial forces and featuring coupling between the bending and torsional vibrations, were investigated. The double beams were connected by a continuous elastic translational and rotational springs. The Chebyshev spectral collocation method was used, and the governing PDEs of motion were discretized. Then, the resulted eigenvalue problem was solved to obtain the critical buckling loads and the natural frequencies. The effects of the axial forces, the elastic layer parameters, the coupling between the bending and torsional vibrations, the slenderness ratio, and the boundary conditions on the critical buckling loads and the natural frequencies were studied. The obtained results show that the bending–torsion coupling has a declining effect on the natural frequencies and the critical buckling loads of single beams, whereas it has a decreasing or increasing effect on the natural frequencies of double beams. Moreover, the veering phenomenon between the frequencies of double beams was observed as a result of varying the axial loads.