Dynamic analysis of axially loaded cantilever shear‐beam under large deflections with small rotations

The modal analysis and seismic response of a vibrating cantilever, with or without tip mass and rotary inertia, are investigated in this study using a shear deformable beam model and including the effect of vertical load. Based on existing approaches, an original method is proposed that does not use fourth‐order uncoupled equations to determine modal deflection and rotation. In fact, the approach presented herein transforms the second‐order coupled system into a first‐order system which can be solved more easily using matrix algebra and Laplace transform. Furthermore, the proposed form allows a straightforward demonstration of orthogonality conditions, that is, the problem is self‐adjoint, and the solution in the case of forced response using modal superposition. In addition, even if the solution presented herein is applicable only to the cantilever with a tip mass and rotary inertia, the scope is general, and the approach can be applied to shear deformable beams with other boundary conditions. Finally, the seismic response by modal superposition is shown, and some examples are proposed and discussed for the case of uniform or continuously varying cross‐sectional properties.


NOVELTY
In this paper, an original approach is used to propose a complete modal analysis for the vibrating compressed cantilever Timoshenko beam with or without tip mass having a rotary moment of inertia, assuming large deflections. Instead of solving the fourth-order uncoupled equations for modal deflection and rotation, in the approach presented in this paper, the second-order coupled system is rewritten in a first-order form which can be solved more easily by resorting to matrix algebra and Laplace transform. A closed-form solution for modal shapes is derived and applications to seismic analysis are shown. the frequencies. A more advanced model was proposed by Rayleigh, which included the effect of kinetic energy due to rotation. However, a significant advancement in beam theories was reflected in the model proposed by Timoshenko 6 , which includes rotary inertia and shear deformability. As highlighted recently, 7 a contribution in the development of the model was given also by Ehrenfest. 8 This model is used in this study. Notably, Traill-Nash and Collar 9 derived the frequencies and modes of a uniform beam without axial force for three types of support. For the same problem, Dolph 10 proposed a solution for the cases of pinned-pinned and free-free end conditions, and established orthogonality conditions. Modal solutions derived differently were proposed by Huang,11 who used the two fourth-order differential equations, that is, one for bending rotation and one for deflection, to obtain frequencies and modes. Abramovich and Elishakoff 12 derived the frequency equation for 10 sets of boundary conditions while considering the individual contributions of shear deformation and rotary inertia but omitting their joint contributions. More recently, Han et al. 5 provided a complete treatment on the evaluation of frequencies and modal shapes based on four beam models, that is, Euler Bernoulli, Rayleigh, Timoshenko and shear-beam (with shear deformability but without rotary inertia). This paper discusses orthogonality so as the response by modal superposition. None of the aforementioned studies include the axial effect, which is relevant in many practical engineering situations.
The effect of the axial force in a Timoshenko beam is included in the study of Kounadis,13 which considers a follower force. Sato 14 derived the complete equations of motion for a cantilever beam in the presence of a tip force directed in the beam-deflected axis. He used the extended Hamilton's principle and assumed that the angle of the deflected axis line with respect to the undeflected line was sufficiently small such that the cosine of that angle was approximately one. This implies that large deflections were considered with small rotation angles; therefore, the difference between the force and its vertical component was negligible.
The same equations were used by Abramovich 15 to derive the characteristic equation for an axially compressed Timoshenko beam under 10 couples of boundary conditions. The solution was obtained by integrating the uncoupled fourth-order equations for deflection and rotation and by imposing the boundary conditions by which deflection and rotation were coupled. No closed-form solutions for the modal shapes were presented.
Another noteworthy problem in practical engineering is a cantilever with a lumped mass and rotary inertia at the free end. The Euler-Bernoulli beam, even in the non-homogeneous case, has been applied in many works, 2,16-19 whereas Timoshenko beams were considered in other studies. [20][21][22] Manolis and Dadoulis 17 proposed a study on the dynamics of a pylon with an attached mass using an Euler-Bernoulli beam with axial deformability. Additionally, they proposed an alternative method to consider the lumped mass directly in the equations of motion, where the Dirac delta function was used instead of applying the mass contribution in the boundary conditions. This approach, although mathematically complicated initially, allows the problem to be solved easily and thus is applied in this study.
This study uses an original approach to propose a complete modal analysis for a vibrating compressed cantilever Timoshenko beam (with or without a tip mass) with a rotary moment of inertia and the assumption of large deflections. The compressive loads are assumed to act vertically, that is, no follower forces are considered. In all the approaches mentioned in the literature, to the best of the author's knowledge, fourth-order uncoupled equations were used to determine the modal deflection and rotation. Subsequently, the constants were determined based on boundary conditions. After fully non-dimensional equations are derived, the approach presented herein transforms the second-order coupled system into a first-order system, which can be solved more easily using matrix algebra and Laplace transform. Instead of being considered in the boundary conditions, the lumped mass is included in the equations using the Dirac delta function, as was performed for the Euler-Bernoulli beam by Manolis and Dadoulis. 17 The proposed form allows a straightforward demonstration of orthogonality conditions, that is, the problem is self-adjoint, and the solution in the case of forced response using modal superposition. In addition, the closed-form expression of modal shapes for uniform beam properties is derived, and F I G U R E 1 Schematic representation of the cantilever beam and reference system (A); internal forces (B). a Galerkin approach is presented for non-uniform beam properties. The solution presented herein is applicable only to a cantilever beam with a tip mass, but the scope is more general, which allows the approach to be applied to beams with other boundary conditions. Finally, the solutions by modal superposition for the seismic base excitation of bridge piers with uniform and variable cross sections are presented and compared with finite element (FE) solutions.

EQUATION OF MOTION
Let us consider a vertical cantilever beam with a reference system (as depicted in Figure 1A), subjected to transversal forces (x,t), vertical time-independent distributed forces (x) (positive in tension based on the reference system in the figure, although we assume that they are always negative), inertia transversal forces, and rotary moments. The equations of motion representing the dynamic equilibrium in the displaced configuration, written in terms of deflection ( , ) and bending rotation ( , ), can be derived as reported by Sato 14 using the extended Hamilton's principle. The fundamental hypothesis is that the angle of rotation of the deflected axis line with respect to the undeformed one is sufficiently small such that the approximations ≅ 1 and ≅ are acceptable, that is, we assume large displacements but small rotations. Notably, is composed of bending rotation and a shear deflection angle. Analogously, the equations can be derived under the same assumptions based on equilibrium and compatibility considerations, 23 as discussed briefly herein.
Concerning the beam element illustrated in Figure 1B, the components of the internal shear and axial forces S(x,t) and N(x,t), respectively, can be related to the components directed along the vertical and transverse directions, P(x) and Q(x,t), respectively, as follows: Considering the previously mentioned hypotheses regarding the deflection angle, we obtain: Considering the equilibrium of the beam element ( Figure 2B) we obtain: where ( ) is the mass, and ( ) = ( ) ( )∕ ( ) is the mass moment of inertia per unit length, ( ) the cross-section area, and ( ) the moment of inertia. We assume that the structure is statically determinate for vertical actions, thus, P(x) is known. Therefore, only the first and third equations of Equation (3) are relevant. The compatibility and constitutive equations lead to the following: Here,̃is the shear rigidity (̃is the area reduced by the shear factor), and is the flexural stiffness. Substituting the first equation of Equation (4) Compared with the equations reported by Sato 14 , Equation (6) contains a vertical internal force term, whereas in the cited paper, the force tangent to the deflected axis is considered; they are equivalent within the assumed hypothesis of a small rotation angle.
The boundary conditions at the free end, in terms of the horizontal component of the internal action and moment, are: The case in which a lumped mass̃and related rotary inertiãare present at the top can be addressed by modifying the boundary conditions as follows: [ However, to verify that the system is self-adjoint such that the equations of motion can be decoupled in the modal space, it is more convenient to consider the mass extremely close to the boundary but internal to the domain of spatial integration. Thus, the boundary conditions shown in Equations (8) In Equations (12) and (13) above, (⋅) is the Dirac delta function that allows us to apply the lumped mass inertia force at = . This is merely a mathematical artifice to achieve our goal. Equations (12) and (13) can be written in nondimensional forms, that is,̂[ using the following normalized quantities: where 0 is the reference frequency, which is assumed to be unitary, and̃( ) is the reference value for shear stiffness. In the case of seismic motion at the base, as will be shown in Section 6, the dimensional force is The boundary conditions can be rewritten using non-dimensional quantities as follows: For completeness, the boundary conditions to be used if Dirac delta is not introduced are as follows: As far as we said, only Equation (18) are used, whereas and are considered in the equations of motion using the Dirac delta function.

MODAL ANALYSIS
Omitting the hat to simplify the non-dimensional variables, the following compact form for the homogeneous problem can be written: where displacement and rotation are compacted in vector form ( , ) = ] and the differential operators  and  are introduced, which are expressed as: Subsequently, solutions are obtained via variable separation as follows: where the same time evolution is assumed for the two modal functions ( ) = [ ( ) Φ( ) ] related to the displacement and bending rotation. Substituting Equation (22) into Equation (20) and denoting the x-derivative by prime for simplicity yields: Based on the first and second equations, the x-and t-dependent functions can be separated as follows: For the previously equations to be valid for all x and t, they must be constant. Therefore, we assumë( )∕ ( ) = − 2 . Thus, the time evolution is governed by:̈( ) + 2 ( ) = 0 (26) with frequency ω determined by imposing boundary conditions. Moreover, the spatial problem is expressed as: The latter is a second-order system of linear differential equations, with variable coefficients, which can be more easily solved after reduction to a first-order system. Finally, previous studies pertaining to shear deformable beams, as cited in the introduction, do not adhere to the procedure shown herein to determine the modal shapes.
Let us define two new variables, ( ) = ′ ( ) and Ψ( ) = Φ ′ ( ). Therefore, Equation (27) can be rewritten as In matrix form: Omitting the argument of the functions for simplicity yields Equation (29) can be solved using the matrix exponential approach only if ( ) commutes with its integral. 24 This is applicable to constant or symmetric matrices but not to general cases. When the uniform beam properties are considered, ′ = ε ′ = ′ = 0 and ε = 1. This case is discussed in this section. Meanwhile, some treatments for non-uniform beams are explained in the next sections.
As will be shown later, it is more convenient to separate the constant component of the matrix from the contributions related to the lumped mass. This can be achieved by expressing as the sum of the following matrices: where: Before addressing the complete problem, let us consider the case without a concentrated mass ( = = 0). Using the matrix exponential, the solution is expressed as follows: where = (0) = [ 1 2 3 4 ] denotes a constant vector. According to the properties of the matrix exponential, the solution can be simplified as follows: In Equation (34), is the matrix whose columns are eigenvectors of matrix 0 , and Λ is the diagonal matrix of eigenvalues, such that 0 = Λ −1 . Thus, the exponential is calculated as for scalar quantities. The eigenvalues of 0 were derived in an ordinary manner via some manipulations, that is, where: For practical purposes, a and b are real. Thus, the eigenvector matrix and its inverse are written as follows: In the previous matrices, i is an imaginary unit. In addition, , , , , and , are dependent on . The complex eigenvalues and eigenvectors form conjugate pairs. Therefore, using the Euler formula in Equation (34) allows one to write the solution based on real quantities using trigonometric expressions, as will be shown later. In the general case involving a top mass Equation (29) can be rewritten by considering Equations (31) and (34) as follows: Performing a Laplace transform on both sides of the equation above and applying the properties of the Dirac delta 25 yields Pre-multiplying by −1 on both sides of the equation above and applying (s) = −1 ( ) yields: By some manipulations, we obtain: Hence, − is diagonal, and the inverse Laplace transform can be calculated as for scalar quantities. By performing an inverse transformation, we obtain In Equation (42), ( ) is the unit-step Heaviside function, defined as 1 for ≥ 1 and 0 otherwise. Finally, This equation can be expressed in terms of the constants (0) only, that is, eliminating (1), as follows: ( The matrix in brackets on the left-hand side is invertible, and its inverse is: − 2 2 − 2 1 . Additionally, = (i,j = 1,2). By some manipulations, the spatial solution can be written in closed form as follows: As shown in Equations (35) and (37), the complex eigenvalues and eigenvectors of 0 form conjugate pairs. Therefore, using the Euler formula in Equation (45) allows one to write the solution as a function of real quantities using trigonometric functions.
The boundary conditions in Equations (18) allows one to write four equations that permit the evaluation of frequencies and constants , as follows: The solutions are those that respect the nil determinant condition In the equation above, the latter is the characteristic equation whose solution provides the structural frequencies, where a, b, r, and s depend on ω. Vector c satisfies Equation (47) for each frequency obtained. One of the non-zero constants is arbitrary. Assuming 3 = 1, 4 can be obtained as follows:

=
( 2 cos − 2 cosh + sin + sinh ) In the equation above, the latter is the solution for 4 using the third equation of Equation (46), and using the fourth equation, of Equation (46), another expression for 4 can be derived, that is, equivalent to the reported one, which is associated with the characteristic Equation (48).

PARAMETRIC ANALYSIS
Before providing an analysis of forced vibrations, the changes in frequencies with the variation in the involved parameters are briefly discussed in this section. When = 1 and = = 0, that is, no vertical load or rotary inertia is considered, the frequency versus β/γ and γ/μ for the first four modes is as shown in Figure 2. Similarly, a plot of the first period, which can be of practical interest, versus β/γ for various γ/μ values is shown in Figure 3.
The effect of the vertical load ( < 1) on the frequency is shown in Figure 4. load. The critical load can be evaluated based on the existing literature for Timoshenko cantilever beams, or equivalently, based on the approach of this study using Equation (33), after setting ω = 0 in Equation (32). By imposing the boundary conditions, the critical value of is = 4∕(4 + 2 ), from which =̃( − 1). Figure 5 shows the modal deflection and rotation for the case presented in Section 7.1, in the presence (a and b) or absence (c and d) of the tip mass. For the same example, to validate the model against other continuous models, Table 1 shows the frequencies obtained using the proposed method when μ = = 0, that is, the case with no tip mass and rotary inertia, with the solution obtained by Abramovich. 15 More comparisons are presented in Section 7 using FE models. Although the FE models are not exact, they are beneficial for addressing aspects not considered by existing continuum models, such as the non-uniformity of the parameters and the presence of tip mass and rotational inertia. Additionally, they are beneficial for seismic analysis.

ORTHOGONALITY CONDITIONS
To obtain an uncoupled system of modal equations, the differential problem must be self-adjoint, as verified using the approach reported by Han et al. 5 Using the differential operators defined in Equation (21) and the total derivative instead of the partial derivative, the problem is self-adjoint if In general, we assume non-uniform beam properties. By performing integration by parts and some manipulations, we obtain The right-hand side of the first equation in Equation (52) can be written as One can easily verify by comparison that the equality is respected, even when non-uniform properties of the beam are considered. The quantities in the brackets are related to the work at the boundaries, the product of the transversal force by displacement, or the bending moment by rotation.
The second equation shown in Equation (53) is immediately verified by considering the properties of the delta function.
Equations (54) and (56) are the orthogonality conditions that allow the equation of motion to be decoupled into linear elementary oscillators.

FORCED RESPONSE BY MODAL SUPERPOSITION
In the case of a forced response, when only time-dependent transverse forces act on the structure, Equation (20) can be written as follows: The response in terms of the superposition of modal shapes is written as follows: By pre-multiplying the above by T i ( ) and performing integration: Subsequently, using the orthogonality conditions in Equations (54) and (56), the modal equations can be expressed as follows: where: When the tip mass is absent, the quantity [ (1)] 2 + [Φ (1)] 2 in the denominator vanishes.
Finally, in the case of seismic motion at the base, considering Equation (17), Equation (62) assumes the following form: When modal damping is considered with ratio , as in the typical case, the term 2 can be added into the previous equation, as will be shown in the examples in Section 7.
Even if the properties of the beam are not constant, Equations (55) and (56) show that if an exact solution is obtained for the frequencies and modal shapes satisfying boundary conditions, then the equations of motion can be decoupled in the modal space. However, in this case, the closed-form solution of the linear system of equations with the variable coefficients in Equation (29) may be difficult to obtain, depending on the complexity of the functions involved. The Galerkin approach can be used to obtain a solution (Appendix A). The first step is to evaluate the approximating functions by assuming constant properties for the beam and fixing the values of the variable coefficients at a specified section. If this section is at the free end, then the approximated modes are evaluated using Equation (50)-(51) also satisfying natural boundary conditions (see Equation (55)). However, numerical analysis shows that this approach is not optimal for the evaluation of base forces. This aspect is discussed in Section 7.3.

Bridge pier with constant properties
As a first example, we consider a bridge pier with a uniform cross-section. The following properties were assumed: The pier was excited at the base via seismic acceleration, which was set as the east component of the earthquake (magnitude Mw 6.5, Figure 6) recorded at the Accumoli station (Italy) on October 30, 2016 by the Italian Strong Motion Network (RAN), 26 during the Central Italy seismic sequence.
The response evaluated using the proposed method was compared with that of a FE model using 60 beam elements, defined using the commercial code Midas Gen (Figure 7). The formulation of the elements was based on the Timoshenko beam theory, where the stiffness effects of tension/compression, shear, and bending deformations were considered. Linear shape functions were used for displacements, whereas quadratic interpolation functions were used for rotation. 27 The FE analysis was based on a direct step-by-step time integration nonlinear analysis, where a damping matrix with a damping ratio of 0.05 was constructed using the modes of the linear model. In the nonlinear case, large displacements were assumed to consider the effect of the vertical load by imposing initial vertical internal element forces to evaluate the geometric stiffness. The frequencies are listed in Table 2.
The results are presented in Figure 8 in terms of the displacement at the top (a), rotation (b), base shear (c), and base moment (d  continuous model. The curves overlapped perfectly. Notably, the FE analysis was nonlinear, whereas the continuous model analysis was linear. This aspect will be further discussed in the following example.

Steel cantilever with pipe section and tip mass
To Two cases were considered: in the first case, the internal vertical force is equal to P = −100 kN, which corresponds to the weight of the tip mass, and in the second case, the tip mass is present but the vertical force is not considered. In this case, the FE model was set up using Midas Gen and 50 beam elements. Two analyses were performed under seismic base excitation using the same time history as that in the previous example. One was based on a linear modal step-by-step integration with 0.05 added modal damping ratio, and the second was a nonlinear analysis, as in the previous example, where large displacements were assumed to consider the effect of the vertical load. The extruded model is shown in Figure 9A. In both the linear and nonlinear cases, a tip mass with rotary inertia as well as shear deformability were considered. Figure 9B-D show the results, as indicated by continuous lines, of the linear FE model (red) and the solution obtained by the proposed method when the vertical force was not considered (black). Moreover, as indicated by dashed lines, the FE nonlinear solution (red) was compared with the proposed method (black) when the force was considered. The solutions for the displacement at the free end, base shear, and base moment are presented. The linear FE case agreed well with the continuous approach solution when vertical force was not considered. By contrast, the nonlinear FE model in large displacements agreed well with the continuous solution with P = −100 kN. Notably, the proposed model is linear, whereas similar FE results are based on a nonlinear model. Table 3 lists the frequencies evaluated by FE linear model and the continuous model.

Bridge pier with variable cross-sections
Finally, a non-uniform cross-section case was considered: a reinforced concrete bridge pier with a hollow cross-section ( Figure 10). The height of the pier was L = 75 m, and the cross-sectional shape, which was almost square, had an external side length that varied based on the following parabolic expression: 9.225 − 0.096 + 0.0006 2 (m). The thickness was 0.6 m. The distributed mass was ( ) = 10 3 (52.832 − 0.587 + 0.004 2 ) (kg), and the vertical variable internal force P = −37751.2 + 518.434 x − 2.88 x 2 + 0.012 x 3 (kN). Figure 11A,B reports the area and moment of inertia versus vertical coordinate. To apply the Galerkin approach described in Section 6 and Appendix A, some approximation functions were evaluated based on the assumption of uniform properties. The cross-sectional properties of the free end were considered in the first analysis. Therefore, the evaluation of in Equation (A.2) was facilitated because the natural boundary conditions were satisfied and the quantities related to the boundaries in Equation (54) vanished.
However, the numerical analysis and comparison with the FE model, based on 75 beam elements showed excellent match for the displacement and rotation at the top, as well as for the base moment, whereas the base shear was underestimated. Therefore, a second analysis was performed, where the functions evaluated based on the properties of the base cross-section were assumed. Only the latter results are reported herein. In this case, was evaluated by considering the boundary conditions in Equation (54). In fact, they were nil at the base, that is, where the essential boundary conditions were applied. By contrast, they did not vanish at the free end when natural conditions were applied. Figure 12 shows the first five modal shapes, that is, the approximated (continuous line) and "corrected" (dashed line) shapes. Table 4 shows the circular frequencies, that is, the approximated circular frequencỹand corrected circular frequency , as well as those derived by FE model. Finally, Figure 13 show a comparison between the FE model (continuous red line) and the continuous model (dashed black line). The agreement shown was reasonable, similarly for the base forces.

CONCLUSIONS
A complete modal analysis was proposed for a compressed cantilever Timoshenko beam with tip mass/inertia, considering only motion in the transverse direction. A continuous model was used for this purpose. An original approach was presented herein for deriving the characteristic equation and modal shapes in a closed form based on a non-dimensional formulation. Instead of solving the fourth-order equations for deriving the modal functions and frequency equation, the proposed approach uses a first-order matrix exponential method. When the tip mass and inertia were considered, the solution was obtained using the Laplace transform. Orthogonality conditions were derived, and the decoupling of the equation of motion in the case of uniform beam properties was presented, such as the solution by modal superposition in the case of a forced seismic response. The Galerkin approach was used when the cross-sectional properties vary with continuity. Comparison with the FE model showed excellent agreement between the response of the continuous linear model under vertical force and that of the nonlinear FE model assuming large displacements. and are non-symmetric n × n matrices, such as −1 . Vectors̃j are not orthogonal with respect to the operators  and . Nevertheless, −1 is diagonalizable and can be written as a product ΔΩ 2 Δ −1 , where Ω 2 is the diagonal matrix