Dynamic Analysis of Timoshenko Beam with Arbitrary Constraints and a Further Optimization Based on Least Energy Principle

Based on Timoshenko’s beam theory, this paper adopts segmented strategy in establishing the governing equations of a multibeam system subjected to various boundary conditions, inwhich free, clamped, hinged, and elastic constraints are considered.Meanwhile, Galerkin method is incorporated as a competitive alternative, in which a new set of unified, efficient, and reliable trial functions are proposed. A further optimization in regard to boundary distributions under forces is implemented and established on the least absorbed energy principle. High agreement is observed between the analytical results and the FEM results, verifying the correctness of the derivations. Complete comparisons between the analytical and the numerical results indicate the Galerkin method is beneficial when slender ratio is larger than 30, in which the continuity of the deformation is proved to be a crucial influencing factor. A modified numerical strategy about optimal boundary is employed and the remarks imply the algorithm can be availably used to reduce the energy absorption of the whole system.


Introduction
The static and dynamic properties of the beam structures in engineering possess great significance, in which the achievements can be directly implemented in wide fields including civil engineering, composite material manufacture, and wing design.In beam analysis, the neglect of the shear and the rotary inertia restricts the utilization of the Euler-Bernoulli theory in short beam cases [1], while the Timoshenko theory can solve these issues well and provide far more reliable solutions [2][3][4].Papers reported in [5,6] were such examples, in which the mode superposition method and the Laplace transform method, for the analytical characteristics, have become the leads.Nevertheless, in some cases that the analytical solution cannot be derived, for example, when system incorporates complicated boundaries, the typical numerical methods like the conventional residual method [7], the Rayleigh-Ritz method [8][9][10], and the Galerkin method [11,12] can provide alternative ways.Further, on account of the great computing ability in computers, the finite element method (FEM) and the boundary element method (BEM) have been developed to calculate more sophisticated structural systems, for example, when system subjected to complicated loads and boundary conditions, which are also regarded as criterions to verify analytical and numerical methods.
The beam systems with various boundary conditions (abbreviated as BCs) and internal constraints (also called BCs), for example, clamped, free, hinged, and elastic BCs, are named as multibeam systems.Actually, such systems have been investigated constantly during recent years [13][14][15].Reference [5] has computed a beam system with hingejoint constraint by employing BEM and FEM and compared the Euler-Bernoulli and the Timoshenko models, in which uniform descriptions for deflection, rotation, bending moment, and shear forces were derived at the first time.Reference [16] built up a FEM model in which the natural frequencies and the mode shapes of an Euler-Bernoulli pontoon system applied with arbitrary numbers of internal hinges were solved via the transfer matrix method.Through a limited superposed method in Euler-Bernoulli model, [17] studied the vibrational power transmission in a uniform and multisupported beam.Reference [18] calculated the eigenvalues by a proposed algorithm and the validity was verified through case studies.Reference [19] considered a multispan Euler-Bernoulli beam with a series of one-sided inner constraints to study the numerical simulations through assumed mode method.Further, [20] analyzed the multispan beams carrying multiple spring-mass systems with axial force based on Timoshenko theory.The free vibration of a multispan Timoshenko beam was studied by Rayleigh-Ritz method in [21], where the trial functions were composed of a series sine or cosine functions plus a set of polynomial functions.Furthermore, by using the transfer matrix method, [22] calculated the eigensolutions of a multispan beam with arbitrary numbers of flexible constraints and compared with the result of the Euler-Bernoulli beam.Reference [23] investigated the free vibrations of the Timoshenko beam with an internal hinge and computed the optimal location of this internal hinge maximized the fundamental frequency.Nevertheless, although most types of the constraints are mentioned in above literatures, few efforts has been made in counting more complex BCs.Additionally, a majority of references lack parallel comparisons with analytical method, the numerical approach, or the FEM.Hence, the primary goal of this article is to analyze a constructed multibeam system established in Timoshenko theory with arbitrary types and numbers of BCs, in which remarkable efforts have been made in derivations and case studies.Based on that, the optimization for BCs under given loads is put forward through the minimum energy principle.This paper is organized as follows.In Section 2, the nondimensional governing equations of the multibeam system are derived through Hamilton's principle and an elementary solution is obtained where several typical BCs are also listed.Section 3 performs the detailed derivations about the beam under arbitrary BCs through an analytical segmented method, in which the matrices are connected through transformation.As an alternative, Galerkin method is introduced in which a group of new trial functions with orthogonal properties are conducted.Section 4 puts forward the optimization process for optimal BCs based on Section 3, achieving the minimum energy absorption under specific loads.A total analysis is drawn in Section 5, in which validity, flexibility, and efficiency of the analytical and the Galerkin methods are measured through FEM in case studies.Based on that, a further optimization about BCs is conducted to verify the derivation in Section 4. Finally, the most relevant conclusions are summarized in Section 6.

Modeling Based on Timoshenko Theory
As depicted in Figure 1, an ordinary beam system is configured under external loads with various BCs (at , the coordinate is   ), each BC herein can be free, simply supported, clamped, hinge-joint, or elastic constraint.As described,   and   denote the external force and moment respectively. is the longitudinal length of the entire beam.Suppose the numbers of elastic BCs are  0 (0 ≤  0 ≤ ) and Hamilton's principle is generally used to construct the governing equations of such multibeam system based on Timoshenko's beam theory.Since this approach is widely used [24][25][26], we refine the derivations as follows.
In consideration of the lateral rotation and the shearing effect, the kinetic energy  and the internal energy  can be expressed as Herein, , , and  represent the bending deflection, bending rotation, and shear deformation respectively, and  and  are the elastic modulus and the shearing modulus. is the socalled shear coefficient [27,28], giving the shear strain ratio on cross section. is the moment inertia in cross section,  is the solid density,   and   are the elastic and torsional spring coefficients of the BC at   , and () is the Dirac delta function.In view of the inherent relationship within Timoshenko model [24], the shear in beam yields For simplification, distributed loading functions in regard to the external forces are originally used and hold In Hamilton's manner,  ∫ 2 1 ( − ) + ∫ 2 1    = 0, and through the variation procedure herein and in conjunction with (2), the final governing differential equations are expressed as For unification, we nondimensionalize (4a) and (4b) by introducing  = /,  = /,  = √/,   = /,  =

Analytical Method (A-Method) through Segmented Strategy.
Considering the big challenges in obtaining the analytical solution of this multiple systems with plenty different BCs, we observe, owning to the segmented geometries, one can solve the vibrational issues via (A.3) for each subdivided part.This inspiration seems falling in between the range of the matrix-transforming method and the FEM [29], but it is generally different from both of them.On the one hand, the arbitrary internal BCs will result in various transfer matrices.On the other hand, it is similar to the solving process of FEM but without any assumptive approximation.Therefore, the strategy to solve this multisystem with  segments proposed here is that (the numbers of BCs are  + 1 and we use  = (0, 1, 2, . . ., ) to express it) we use the fundamental solution shown in Appendix A to solve   = ( ) first and connect them in sequence via physical boundary relations, by which we are able to get the nature frequency   and its modal shape Φ  through eigenanalysis.Herein, a classified discussion about different BCs is listed as follows due to the arbitrariness of the BCs.
At first, a series of dimensionless transformations in localised coordinates are made here.For segment , we introduce   = (  −  −1 )  , X  = ( −  −1 )/(  −  −1 ), Z  = /(  −  −1 ), and   =  to guarantee a uniform  throughout the segments, and as a ripple effect,   and   in each fundamental solutions are changed into Hence, the shape function of segment  naturally holds In which    = (cosh(    ), sinh(    ), cos(    ), sin(    )).To get the specific Φ  (X  ), each   is to be determined.In view of the dissimilarity of BCs at   , we generate the following: (1) When  = 0 (BC at  = 0), we enumerate the decompositions of the BCs as follows: (a) For free BC Mathematical Problems in Engineering (b) For simply supported BC (2) When 0 <  < , the free constraint disappears, but the hinge-joint BC is added.
(a) For simply supported BC (3) When  = , accordingly one has the following: (a) For free BC In conclusion, all mentioned expressions have generated two linear equations at two edge nodes and four linear equations at inner nodes about C (detailed information is shown in Appendix B), which constitute 4 ×  linear equations of  nodes in total denoted by where the coefficient matrix B just involves the independent unknown variable .We further define the characteristic function as By which we can obtain each   in sequence after solving this super sophisticated transcendental equation Δ 1 = 0.Meanwhile, C  and Φ   (X  ) are available to obtain through Ultimately, for  −1 ≤  ≤   , through coordinate transformation, we generate

Numerical Solution Based on Galerkin Method (G-Method).
Albeit the achievement in analytical method, the explicit Δ 1 is far more complex and costly in actual computation, especially when  becomes larger.Actually, in allusion to the practicability, the static and dynamic properties of the multisystem, in many circumstances, are determined by the lower modes ( = 1, 2).Based on such consideration, the numerical solution, as an alternative way, may probably provide credible results with avoiding massive computing costs.Hence, the objective of this part is the derivations about numerical G-method.
Based on G-method, the residual functions of the entire system are where L and B are all differential operators about the physical residual and the natural residual. * and  * separately, as the external forces and the physical constraints of BC, stand for the given restricting equations independent on the unknown function Φ. Herein,  * = 0 is naturally satisfied.In general, the trail functions possess Accordingly, the functions of residual errors are To reduce the complexity,   = 0 is usually satisfied in advance [11,22,25], which implies that all G 1 1 () and G 1 2 () should be constructed at first.According to Appendix C, the explicit G 1 1 and G 1 2 for free BCs are hard to work out because of the coupling in (C.1) about displacement and rotation from shear component.Nevertheless, if we reduce the demands, that is, to decouple (C.1) and replace it by a weaker Euler-Bernoulli condition in (27), not only will the difficulty be solved, but also the inner beam maintains Timoshenko precision To guarantee an entire   = 0, the polynomials are employed to structure the fundamental trial functions and a series of G 1 1 and G 1 2 are proposed by expanding previous result respectively via polynomial, Legendre orthogonal polynomials, Chebyshenv orthogonal polynomials, and trigonometric functions, which is regarded as the explorations of preferable Φ trial .Here, we have given the exact nature conditions in Appendix C that the trial functions should hold.Further, if the system has not clamped BCs in 0 <  <  (actually, if clamped BCs exist in inner nodes, the total system can be divided into several individual beams), a set of explicit trial functions can be proposed as which include four kinds of F  : (1) Polynomial form (P-method) (2) Chebyshev form (C-method) (3) Legendre form (L-method) With the explicit Φ trial (), the final derivation is obtainable, and we conduct and make After substituting the specific expressions to (34) and a group of integrations [11,12] (detailed in Appendix D), the special variables are vanished, leaving behind a time domain equation denoted as The characteristic equation can be easily carried out through Δ 2 = det(− 2 +), by which it is capable of getting the approximate frequencies through Δ 2 = 0.According to ), the corresponding eigenvector can be solved and consequently described as ) . (36)

The Optimization of BCs
Under external loads, the positions of BCs will determine the dynamic performance of the whole system, wherein a reasonable distribution about these nodes may improve the loading capacity and its stable margin [23].Based on the work in Section 3, it is capable of acquiring the optimal BCs.For mode , through a series of modal transformations as reported in [30], we have Under certain loads, we write Φ  (, ) = ( () )   (), in modal superposition manner; (37a) and (37b) satisfy Based on modal orthogonality, we obtain To examine the quantitative performance, we further define the total absorbed energy of the whole system as Hence, to achieve the minimum energy absorption, (40) must yield Θ = Θ min .In virtue of the derivations in Section 3, we get  *  =  *  (,  1 ,  2 , . . .,  −1 ) and   =   ( 1 ,  2 , . . .,  −1 ), all of which are the functions of ( 1 ,  2 , . . .,  −1 ).Herein, we suppose the beam without initial disturbance and, after Duhamel integration, get the modal response as After substituting it to (40), we obtain Particularly for time-invariant loads, it is relatively easy to calculate the absorbed energy when system reaches into equilibrium state ( = 0), which holds Further, to derive the optimal layout, (42) and (43) must satisfy (44), which comes down to a minimization issue By which the optimal ( 1 ,  2 , . . .,  −1 ) can be worked out.

Numerical Verification
To demonstrate the validity of the proposed solutions and make comparisons each other, we conduct 2 typical numerical examples illustrated in Figure 2 named as case 1 and case 2, respectively, in which different BCs are investigated and inherent variables are concerned, for example, the slender ratio   and the BC positions   .For the purpose of evaluation, the results of FEM are generated through Abaqus 6.13, whereby the verification is demonstrated and complete comparisons are made.Further, based on the effectiveness of Section 5.2, case 2 is chosen to carry out a further optimization study in Section 5.3.

Considerations during Computations.
Before the numerations, several considerations are declared to obtain more reliable consequences.(1) Strategy of algorithm: since the numerations in such issues involve massive matrix computations, Matlab programs are utilized.Whereas in special situations, for example,  1 = 0.01 or  2 = 0.02, the bad condition-numbers of  or  may exist and distort the precision.Based on that, an examination and repairing step is incorporated during the calculation.Additionally, we reserve the real part of det(B) to solve some big eigenmatrices and acquire the eigenvalues and eigenvectors, holding real(det (B)) = 0 which is proved with pretty effectiveness.
(2) "Mirror verification" (interpreted as algorithm flexibility analysis): for a set of symmetrical BCs, for example,  1 = 0.3,  2 = 0.6 and  1 = 0.4,  2 = 0.7, the results of mode frequencies and mode shapes in each pair should be theoretically identical and symmetrical.This principle is regarded as a measurement about algorithm stability for different methods, whereby we call the "mirror verification." (3) Gridding scheme in FEM: basically, the accuracy of the FEM mostly depends on the mesh and the grid number, and hence the distributions of meshes as depicted in Figure 2 guarantee that no matter how long the proportion in each part, as sufficient as 10000 nodes are generated for each segment to improve the reliability of FEM.

Cases Studies.
As shown in Figure 2, a beam with two elastic edges and two internal hinged joints are considered in case 1, while in case 2 a beam with two free BCs and two simply supported internal constrains are investigated.Apparently, the rotational deformation in case 1 shows discontinuity, while case 2 is continuous as the opposite situation.Based on above specific configurations, it is able to construct the definite governing matrices B according to (8a) to (19c) and (28a) and (28b).Herein,  = 3 is employed in the study.
And for G-method, it holds (1) case 1  (2) case 2 During actual computations, ( 21) and ( 23) are applied to get modal frequencies and modal shapes.To solve the big 12 × 12 eigenmatrix of ( 22), the proposed algorithm strategy is applied.Meanwhile, the gridding scheme is also applied to the FEM.
As presented in Figure 3, the frequency plots of FEM and A-method at 1 = 0.3 and 2 = 0.6 are depicted, respectively, in which three slender ratios are incorporated as the representatives of general short to long beams [1].Results show that, no matter in case 1 or case 2, A-method shares same tendencies with FEM in which high coherence is observed in various   .Figure 4 gives the corresponding modal shapes through normalization in each dimension, where only the natural displacement components are illustrated.It is observed that perfect coincidence appears except for the 1st and 3rd modes in case 1, demonstrating the overall correctness of the A-method.
Although the comparisons with FEM obtain satisfactory results, the high costs of A-method in computation, where the computation under  = 3 even occupies tremendous costs, limits its application.Therefore, a total comparison is made between the G-method and the A-method to explore the potential substitutability.Figure 5 displays the outcome of the four G-methods based on polynomial, Chebyshev, Lagrange, and trigonometric functions, respectively, in which the computations consider a set of symmetric BC positions ( 1 = 0.3,  2 = 0.6 and  1 = 0.4,  2 = 0.7) and various  1 and  2 are utilized.The comparisons within case 1, as exhibited in Figures 5(a) and 5(b), show there exist apparent errors between the A-method and the Gmethod.Wherein the 1st frequency generated by G-method lies between the 2nd and the 3rd modes calculated through A-method, which can be interpreted as, although the trial functions used in (13a), (13b), (13c), and (13d) already hold   = 0, the modal sectional angles at hinge-joint BCs ( 1 and  2 ) as demonstrated in Figure 6(b) inherit discontinuities, which exceeds the general ability of (13a), (13b), (13c), and (13d) and results in imprecise results.In case 2, the results turn better, as the displacement  and the rotation  at fixhinge BCs transform smoothly (as shown in Figures 6(c) and 6(d)), and the roots of 1st frequency except for T-method fall in the adjacent range of exact value uniformly as pictured in Figures 5(c) and 5(d).The quantitative computing effect of G-method can be further measured through correlation coefficients relative to A-method [31].
Besides, it is undeniable that the slender ratio is regarded as an essential parameter in beam system [1,32].Hence, a large-scale variations about   from 1 to 100 are carried out.Overall, we capture in Figure 7 that the curves of frequencies present irregular changes due to the discontinuity in case 1, while the trends are relatively uniform in case 2. Detailed information observed in Figures 7(a) and 7(b) shows that the numerical methods show surging changes in the plots within   < 25, demonstrating the numerical methods is not absolutely perfect for all slenderness and there still exist certain algorithm instabilities.For case 2, as illustrated in Figure 7(d), the C-method shows break-offs when   > 26 which indicates the P-method and the L-method are   more stable than the C-method and the T-method.Actually, the entire curve tendencies in Figure 7 demonstrate that the numerical precisions are significantly dependent on   , where bad computing effect occurs within low   , while in large   the consequences of G-methods are fairly acceptable.Further, the curves inherit closing regions near   = 15 to   = 25 in case 1, after which the solving effect of G-methods slightly decreases.In case 2, although the G-method and the A-method share the consistent trends, the solving effect is indeed affected by   .On the whole, we conclude the Pmethod and the L-method possess higher flexibility.It is also indicated that the G-methods can be deemed with satisfied precision compared with the A-method when   > 30.

BCs' Optimization Based on Case 2.
Since the numerical method can predict the general frequency tendencies with reasonable and low-consumed properties, a further optimization based on case 2 is carried out in this part.Actually, the objective of the optimization conducted here is to obtain the optimal BCs proposed in Section 4 to achieve the minimum energy reservation in actual structural design.The A-method and the G-method are also applied together as a contrast.It is undeniable that the explicit expression of ( 46) is hard to get even for case 2 with only 2 fixed-hinge BCs.On the one hand, although the det step performed in (21) can obtain the characteristic equations with variables  1 and  2 , to solve the frequency is almost impossible with two unknown parameters.On the other hand, the derivations of (46) about  1 and  2 have introduced a set of implicit characteristic equations, which are considered to be super complicated.Thus, a numerical strategy is put forward here on the base of Section 4.
Algorithm scheme: as shown in Figure 8, we traverse  1 and  2 in as small as possible step simultaneously to get reliable outcome during numerations and form the initial database.On such basis, a fitting procedure is drawn through Matlab's griddata(. . ., "V4" ) function (herein, "V4" means 4 grid points spline functions interpolation).In this process, 1000 × 1000 mesh grids (1 million combinations in total) are generated for  1 and  2 , after which the minimum energy interpolating point is sought out.Sequentially, the dualistic algebraic interpolations are expanded at this point within the chosen grid area denoted as Θ * [33].Afterwards, a procedure of derivations in (46) about  1 and  2 is made and the twodimension equations are conducted to get the final optimal BC values.Actually, excessive integrals must be solved at each  1 and  2 combination during the traversal process to get the specific energy, which is considered to be a computationally expensive procedure.The balance of the traversal step must be considered due to the massive computing time, especially for A-method.
To conduct the numeration, we set   = 10 −6 / and   = 0. To evaluate the numerical methods, the reservation of the first 20 modal energies in A-method is regarded as the  accurate solution.According to the computing performance, we choose 0.05 in A-method, 0.01 in P-method, 0.05 in Cmethod, 0.02 in L-method, and 0.02 in T-method as the inherent steps for  1 or  2 , respectively.Figures 9 and 10 are the frequency spectrums and the energy maps of the Amethod and the G-methods, in which the frequencies and the energy clouds are in keeping with each other.It is found that the BCs at higher frequencies tend to absorb more external energies generally in which, although the peaks and troughs are irregular, the G-methods perform good accordance to a certain degree.On the whole, we observe from Figures 9 and 10 that the P-method and the L-method present good "Mirror" properties, and the C-method is slightly inferior, while the T-method performs worse.The final optimal BCs of the two methods are listed in Table 1 in which the errors indicate the G-methods can effectively predict the optimal positions, demonstrating the substitutability of the numerical method.

Conclusion
Based on Timoshenko's beam theory, this paper studies the vibrational performance of a continuous beam subjected to     arbitrary types and numbers of constrains.Based on that, the optimization about boundary conditions, established on the minimum energy absorbed principle under loading, is proposed.During the derivations, analytical method is employed at first based on segmented strategy and the numerical Galerkin method is also adopted for comparison and as an alternative approach, in which a new type of trial functions is constructed.Further, the analytical equations are derived for boundary optimization considering the external loads.In case studies, two typical numerical examples representing the discontinuous and the continuous situations, respectively, are introduced.The results indicate that the analytical method is valid compared with the FEM but costly indeed.The comparisons between the analytical method and the Galerkin methods demonstrate the precisions of the numerical method are higher in continuous case, since the discontinuity in mode shapes hinders the entire simulation accuracy.It is also verified the Galerkin method is beneficial in the case of   > 30.Furthermore, a practical optimization is conducted based on case 2, in which a suitable computing strategy is put forward and the effectiveness of the G-methods are compared with

C. Boundary Conditions in Trial Functions
The nature conditions of trial functions in Galerkin method are as follows.For free boundaries ( = 0 or  = 1), the trail functions must yield

Figure 4 :
Figure 4: The first 12 displacemental modal shapes of A-method and FEM.

Figure 5 :
Figure 5: Frequencies generated by different G-methods with various  1 and  2 .
angles between A-method and Galerkin (d) Case 2

Figure 7 :
Figure 7: Frequency curves of A-method and G-methods with change of   .

Figure 9 :
Figure 9: Frequency spectrums computed through A-method and G-methods.

Table 1 :
Optimizing BCs and errors compared with A-method.