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

Authors: Xiayang Zhang; Ming Zhu; Zhe Wu; Meijuan Zhao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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, in which 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.


---

## Body

## 1. 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–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–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–15]. Reference [5] has computed a beam system with hinge-joint 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 Section2, 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.


## 2. Modeling Based on Timoshenko Theory

As depicted in Figure1, an ordinary beam system is configured under external loads with various BCs (at i, the coordinate is xi), each BC herein can be free, simply supported, clamped, hinge-joint, or elastic constraint. As described, Fp and Fm denote the external force and moment respectively. L is the longitudinal length of the entire beam. Suppose the numbers of elastic BCs are J0 (0≤J0≤N) 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–26], we refine the derivations as follows.Figure 1
Configuration of the beam.In consideration of the lateral rotation and the shearing effect, the kinetic energyT and the internal energy V can be expressed as(1a)T=12∫0LρA∂zx,t∂t2dx+12∫0LρI∂φx,t∂t2dx,(1b)V=12∫0LEI∂φx,t∂x2dx+12∫0LkGAγx,t2dx+12∑J=1J0ksJzx,tδx-xJ2+12∑J=1J0ktJφx,tδx-xJ2.Herein, z, φ, and γ represent the bending deflection, bending rotation, and shear deformation respectively, and E and G are the elastic modulus and the shearing modulus. k is the so-called shear coefficient [27, 28], giving the shear strain ratio on cross section. I is the moment inertia in cross section, ρ is the solid density, ksJ and ktJ are the elastic and torsional spring coefficients of the BC at xJ, and δ() is the Dirac delta function. In view of the inherent relationship within Timoshenko model [24], the shear in beam yields(2)γx,t=∂zx,t∂t-φx,t.For simplification, distributed loading functions in regard to the external forces are originally used and hold(3)δW=∫0LFpx,tδzx,tdx+∫0LFmx,tδφx,tdx.In Hamilton’s manner,δ∫t1t2T-Vdt+∫t1t2δWncdt=0, and through the variation procedure herein and in conjunction with (2), the final governing differential equations are expressed as(4a)-∂∂xkGA∂zx,t∂x-φx,t+ρA∂2zx,t∂t2+∑J=1J0ksJzx,tδx-xJ=Fpx,t,(4b)-∂∂xEI∂φx,t∂x-kGA∂zx,t∂x-φx,t+ρI∂2φx,t∂t2+∑J=1J0ktJφx,tδx-xJ=Fmx,t.For unification, we nondimensionalize (4a) and (4b) by introducing X=x/L, Z=z/L, λ=kG/E, Lr=L/r, t=τ/w0, w0=kGE/ρL2, Fp∗X,τ=Fpx=XL,t=τ/w0/ρALw02, Fm∗(X,τ)=Fmx=XL,t=τ/w0/ρIw02, ksJ∗=ksJ/ρALw02, and ktJ∗=ktJ/ρILw02. Thus, (4a) and (4b) are rewritten as(5a)Z¨-λZ″+λφ′+∑J=1J0ksJ∗ZδX-XJ=Fp∗,(5b)φ¨-1λφ″+λLr2φ-λLr2Z′+∑J=1J0ktJ∗φδX-XJ=Fm∗,where symbols   ¨,   ˙, ″, and ′ are employed to denote ∂2/∂τ2, ∂/∂τ, ∂2/∂X2, and ∂/∂X, respectively.


## 3. Derivations of Solution

### 3.1. 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 N segments proposed here is that (the numbers of BCs are N+1 and we use i=(0,1,2,…,N) to express it) we use the fundamental solution shown in Appendix A to solve Ci=Ci1Ci2Ci3Ci4ki2Ci2ki1Ci1ki4Ci4ki3Ci3 first and connect them in sequence via physical boundary relations, by which we are able to get the nature frequency wj and its modal shape Φj 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 segmenti, we introduce Lri=Xi-Xi-1Lr, Xi=X-Xi-1/Xi-Xi-1, Zi=Z/(Xi-Xi-1), and λi=λ to guarantee a uniform w throughout the segments, and as a ripple effect, αi and βi in each fundamental solutions are changed into(6a)αi=12Lri-1-1λ2+1-1λ22+4Lri2λw2Xi-Xi-121/21/2,(6b)βi=12Lri1+1λ2+1-1λ22+4Lri2λw2Xi-Xi-121/21/2.Hence, the shape function of segmenti naturally holds(7)ΦiXi=CiPi.In which PiT=(cosh⁡biαiX,sinh⁡biαiX,cos⁡biβiX, sin⁡biβiX). To get the specific ΦiXi, each Ci is to be determined. In view of the dissimilarity of BCs at Xi, we generate the following: (1)
Wheni=0 (BC at X=0), we enumerate the decompositions of the BCs as follows:
(a)
For free BC(8a)Φ1′2,1X1=0=0,Φ1′1,1X1=0-Φ12,1X1=0=0,(8b)B0I10P1X1=0=0,(b)
For simply supported BC(9a)Φ11,1X1=0=0,Φ1′2,1X1=0=0,(9b)B0I20P1X1=0=0.(c)
For clamped BC(10a)Φ11,1X1=0=0,Φ12,1X1=0=0,(10b)B0I30P1X1=0=0.(d)
For elastic BC(11a)λΦ1′1,1-Φ12,1X1=0=ks0∗Φ11,1X1=0X1,1λX1Φ12,1′X1=0=kt0∗Φ12,1X1=0,(11b)B0I40P1X1=0=0.(2)
When0<i<N, the free constraint disappears, but the hinge-joint BC is added.
(a)
For simply supported BC(12a)Φi+11,1Xi+1=0=0,Φi1,1Xi=1=0,(12b)Φi+12,1Xi+1=0=Φi2,1Xi=1,(12c)Φi+1′2,1Xi+1=0Xi+1-Xi=Φi′2,1Xi=1Xi-Xi-1,(12d)Bi1PiXi=1Pi+1Xi+1=0=0.(b)
For hinge-joint BC(13a)Φi1,1Xi=1Xi-Xi-1=Φi+11,1Xi+1=0Xi+1-Xi,(13b)Φi′2,1Xi=1=0,Φi+1′2,1Xi+1=0=0,(13c)Φi′1,1-Φi2,1Xi=1=Φi+1′1,1-Φi+12,1Xi+1=0,(13d)Bi2PiXi=1Pi+1Xi+1=0=0.(c)
For clamped BC(14a)Φi+11,1Xi+1=0=0,Φi+12,1Xi+1=0=0,(14b)Φi1,1Xi=1=0,Φi2,1Xi=1=0,(14c)Bi3PiXi=1Pi+1Xi+1=0=0.(d)
For elastic BC(15a)Φi1,1Xi=1Xi-Xi-1=Φi+11,1Xi+1=0Xi+1-Xi,(15b)Φi2,1Xi=1=Φi+12,1Xi+1=0,(15c)λΦi+1′1,1-Φi+12,1Xi+1=0-Φi′1,1-Φi2,1Xi=1=ksi∗Φi+11,1Xi+1=0Xi+1-Xiorksi∗Φi1,1Xi=1Xi-Xi-1,(15d)1λΦi+12,1′Xi+1=0Xi+1-Xi-Φi2,1′Xi=1Xi-Xi-1=kti∗Φi+12,1Xi+1=0orkti∗Φi2,1Xi=1,(15e)Bi4PiXi=1Pi+1Xi+1=0=0.(3)
Wheni=N, accordingly one has the following:
(a)
For free BC(16a)ΦN′2,1XN=1=0,ΦN′1,1XN=1-ΦN2,1XN=1=0,(16b)BNI11PNXN=1=0(b)
For simply supported BC(17a)ΦNN,1XN=1=0,ΦN′2,1XN=1=0,(17b)BNI21PNXN=1=0.(c)
For clamped BC(18a)ΦNN,1XN=1=0,ΦN2,1XN=0=0,(18b)BNI31PNXN=1=0.(d)
For elastic BC(19a)λΦN′1,1-ΦN2,1XN=1=ksN∗ΦN1,1X1=1XN-XN-1,(19b)1λXN-XN-1ΦN2,1′XN=1=ktN∗ΦN2,1XN=1,(19c)BNI41PNXN=1=0.In conclusion, all mentioned expressions have generated two linear equations at two edge nodes and four linear equations at inner nodes aboutC (detailed information is shown in Appendix B), which constitute 4×N linear equations of N nodes in total denoted by CT=C1,1,C1,2,C1,3,C1,4,…,CN1,CN,2,CN,3,CN,4. That is,(20)BC=04N×1,where the coefficient matrix B just involves the independent unknown variable w. We further define the characteristic function as(21)Δ1=detB.By which we can obtain each wj in sequence after solving this super sophisticated transcendental equation Δ1=0. Meanwhile, Cj and ΦijXi are available to obtain through (22)Bw=wjCj=0.Ultimately, forXi-1≤X≤Xi, through coordinate transformation, we generate(23)ΦjX=LriLrΦijXi=X-Xi-1/Xi-Xi-1=Xi-Xi-1ΦijXi=X-Xi-1/Xi-Xi-1.


### 3.2. 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 N 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 j=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(24a)L~Φ-F∗=0,X≠Xi,(24b)B~Φ-g∗=0,X=Xi,where L~ and B~ are all differential operators about the physical residual and the natural residual. F∗ and g∗ separately, as the external forces and the physical constraints of BC, stand for the given restricting equations independent on the unknown function Φ. Herein, g∗=0 is naturally satisfied. In general, the trail functions possess(25)ΦtrialX,τ=∑i1=1N1G1i1XQ1i1τ∑i2=1N2G2i2XQ2i2τ.Accordingly, the functions of residual errors are(26a)RI=L~Φtrial-F∗,X≠Xi,(26b)RB=B~Φtrial,X=Xi.To reduce the complexity,RB=0 is usually satisfied in advance [11, 22, 25], which implies that all G1i1X and G1i2X should be constructed at first. According to Appendix C, the explicit G1i1 and G1i2 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(27)G2i2I10X=0″=G2i2I11X=1″=0.To guarantee an entireRB=0, the polynomials are employed to structure the fundamental trial functions and a series of G1i1 and G1i2 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<i<N (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(28a)G1i1=∏I20ifexistX∏I21ifexistX-1∏j1=1ifexist,J1X-Xj1∏j3=1ifexist,J3X-Xj3Fi1X,(28b)G2i2=∫∏I10ifexistX2∏I11ifexistX-12∏I20ifexistX∏I21ifexistX-1∏j2=1ifexist,J2X-Xj2Fi2XdX+Constant,which include four kinds of Fi: (1)
Polynomial form (P-method)(29)FjPX=Xj-1,j≥1,(2)
Chebyshev form (C-method)(30a)F1C=1,(30b)F2C=2X-1,(30c)Fj+1C=22X-1FjC-Fj-1C,j≥2.(3)
Legendre form (L-method)(31a)F1L=1,(31b)FjL=14jj!djdXj2X-12-1j,j≥2.(4)
Trigonometric form (T-method)(32)FjTX=cos⁡j-1X+1,j≥1.With the explicitΦtrialX, the final derivation is obtainable, and we conduct(33)GT=G11G12⋯G1N100⋯000⋯0G21G22⋯G2N1.and make(34)∫01GRIdX=0After 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(35)MQ¨+KQ=F.The characteristic equation can be easily carried out throughΔ2=det⁡-w2M+K, by which it is capable of getting the approximate frequencies through Δ2=0. According to -w¯j2M+KQj=0 (here Q=Qjqτ), the corresponding eigenvector can be solved and consequently described as(36)Φ¯j=∑i1=1N1G1i1XQji1,1∑i2=1N2G2i2XQji2,1.


## 3.1. 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 N segments proposed here is that (the numbers of BCs are N+1 and we use i=(0,1,2,…,N) to express it) we use the fundamental solution shown in Appendix A to solve Ci=Ci1Ci2Ci3Ci4ki2Ci2ki1Ci1ki4Ci4ki3Ci3 first and connect them in sequence via physical boundary relations, by which we are able to get the nature frequency wj and its modal shape Φj 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 segmenti, we introduce Lri=Xi-Xi-1Lr, Xi=X-Xi-1/Xi-Xi-1, Zi=Z/(Xi-Xi-1), and λi=λ to guarantee a uniform w throughout the segments, and as a ripple effect, αi and βi in each fundamental solutions are changed into(6a)αi=12Lri-1-1λ2+1-1λ22+4Lri2λw2Xi-Xi-121/21/2,(6b)βi=12Lri1+1λ2+1-1λ22+4Lri2λw2Xi-Xi-121/21/2.Hence, the shape function of segmenti naturally holds(7)ΦiXi=CiPi.In which PiT=(cosh⁡biαiX,sinh⁡biαiX,cos⁡biβiX, sin⁡biβiX). To get the specific ΦiXi, each Ci is to be determined. In view of the dissimilarity of BCs at Xi, we generate the following: (1)
Wheni=0 (BC at X=0), we enumerate the decompositions of the BCs as follows:
(a)
For free BC(8a)Φ1′2,1X1=0=0,Φ1′1,1X1=0-Φ12,1X1=0=0,(8b)B0I10P1X1=0=0,(b)
For simply supported BC(9a)Φ11,1X1=0=0,Φ1′2,1X1=0=0,(9b)B0I20P1X1=0=0.(c)
For clamped BC(10a)Φ11,1X1=0=0,Φ12,1X1=0=0,(10b)B0I30P1X1=0=0.(d)
For elastic BC(11a)λΦ1′1,1-Φ12,1X1=0=ks0∗Φ11,1X1=0X1,1λX1Φ12,1′X1=0=kt0∗Φ12,1X1=0,(11b)B0I40P1X1=0=0.(2)
When0<i<N, the free constraint disappears, but the hinge-joint BC is added.
(a)
For simply supported BC(12a)Φi+11,1Xi+1=0=0,Φi1,1Xi=1=0,(12b)Φi+12,1Xi+1=0=Φi2,1Xi=1,(12c)Φi+1′2,1Xi+1=0Xi+1-Xi=Φi′2,1Xi=1Xi-Xi-1,(12d)Bi1PiXi=1Pi+1Xi+1=0=0.(b)
For hinge-joint BC(13a)Φi1,1Xi=1Xi-Xi-1=Φi+11,1Xi+1=0Xi+1-Xi,(13b)Φi′2,1Xi=1=0,Φi+1′2,1Xi+1=0=0,(13c)Φi′1,1-Φi2,1Xi=1=Φi+1′1,1-Φi+12,1Xi+1=0,(13d)Bi2PiXi=1Pi+1Xi+1=0=0.(c)
For clamped BC(14a)Φi+11,1Xi+1=0=0,Φi+12,1Xi+1=0=0,(14b)Φi1,1Xi=1=0,Φi2,1Xi=1=0,(14c)Bi3PiXi=1Pi+1Xi+1=0=0.(d)
For elastic BC(15a)Φi1,1Xi=1Xi-Xi-1=Φi+11,1Xi+1=0Xi+1-Xi,(15b)Φi2,1Xi=1=Φi+12,1Xi+1=0,(15c)λΦi+1′1,1-Φi+12,1Xi+1=0-Φi′1,1-Φi2,1Xi=1=ksi∗Φi+11,1Xi+1=0Xi+1-Xiorksi∗Φi1,1Xi=1Xi-Xi-1,(15d)1λΦi+12,1′Xi+1=0Xi+1-Xi-Φi2,1′Xi=1Xi-Xi-1=kti∗Φi+12,1Xi+1=0orkti∗Φi2,1Xi=1,(15e)Bi4PiXi=1Pi+1Xi+1=0=0.(3)
Wheni=N, accordingly one has the following:
(a)
For free BC(16a)ΦN′2,1XN=1=0,ΦN′1,1XN=1-ΦN2,1XN=1=0,(16b)BNI11PNXN=1=0(b)
For simply supported BC(17a)ΦNN,1XN=1=0,ΦN′2,1XN=1=0,(17b)BNI21PNXN=1=0.(c)
For clamped BC(18a)ΦNN,1XN=1=0,ΦN2,1XN=0=0,(18b)BNI31PNXN=1=0.(d)
For elastic BC(19a)λΦN′1,1-ΦN2,1XN=1=ksN∗ΦN1,1X1=1XN-XN-1,(19b)1λXN-XN-1ΦN2,1′XN=1=ktN∗ΦN2,1XN=1,(19c)BNI41PNXN=1=0.In conclusion, all mentioned expressions have generated two linear equations at two edge nodes and four linear equations at inner nodes aboutC (detailed information is shown in Appendix B), which constitute 4×N linear equations of N nodes in total denoted by CT=C1,1,C1,2,C1,3,C1,4,…,CN1,CN,2,CN,3,CN,4. That is,(20)BC=04N×1,where the coefficient matrix B just involves the independent unknown variable w. We further define the characteristic function as(21)Δ1=detB.By which we can obtain each wj in sequence after solving this super sophisticated transcendental equation Δ1=0. Meanwhile, Cj and ΦijXi are available to obtain through (22)Bw=wjCj=0.Ultimately, forXi-1≤X≤Xi, through coordinate transformation, we generate(23)ΦjX=LriLrΦijXi=X-Xi-1/Xi-Xi-1=Xi-Xi-1ΦijXi=X-Xi-1/Xi-Xi-1.


## 3.2. 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 N 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 j=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(24a)L~Φ-F∗=0,X≠Xi,(24b)B~Φ-g∗=0,X=Xi,where L~ and B~ are all differential operators about the physical residual and the natural residual. F∗ and g∗ separately, as the external forces and the physical constraints of BC, stand for the given restricting equations independent on the unknown function Φ. Herein, g∗=0 is naturally satisfied. In general, the trail functions possess(25)ΦtrialX,τ=∑i1=1N1G1i1XQ1i1τ∑i2=1N2G2i2XQ2i2τ.Accordingly, the functions of residual errors are(26a)RI=L~Φtrial-F∗,X≠Xi,(26b)RB=B~Φtrial,X=Xi.To reduce the complexity,RB=0 is usually satisfied in advance [11, 22, 25], which implies that all G1i1X and G1i2X should be constructed at first. According to Appendix C, the explicit G1i1 and G1i2 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(27)G2i2I10X=0″=G2i2I11X=1″=0.To guarantee an entireRB=0, the polynomials are employed to structure the fundamental trial functions and a series of G1i1 and G1i2 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<i<N (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(28a)G1i1=∏I20ifexistX∏I21ifexistX-1∏j1=1ifexist,J1X-Xj1∏j3=1ifexist,J3X-Xj3Fi1X,(28b)G2i2=∫∏I10ifexistX2∏I11ifexistX-12∏I20ifexistX∏I21ifexistX-1∏j2=1ifexist,J2X-Xj2Fi2XdX+Constant,which include four kinds of Fi: (1)
Polynomial form (P-method)(29)FjPX=Xj-1,j≥1,(2)
Chebyshev form (C-method)(30a)F1C=1,(30b)F2C=2X-1,(30c)Fj+1C=22X-1FjC-Fj-1C,j≥2.(3)
Legendre form (L-method)(31a)F1L=1,(31b)FjL=14jj!djdXj2X-12-1j,j≥2.(4)
Trigonometric form (T-method)(32)FjTX=cos⁡j-1X+1,j≥1.With the explicitΦtrialX, the final derivation is obtainable, and we conduct(33)GT=G11G12⋯G1N100⋯000⋯0G21G22⋯G2N1.and make(34)∫01GRIdX=0After 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(35)MQ¨+KQ=F.The characteristic equation can be easily carried out throughΔ2=det⁡-w2M+K, by which it is capable of getting the approximate frequencies through Δ2=0. According to -w¯j2M+KQj=0 (here Q=Qjqτ), the corresponding eigenvector can be solved and consequently described as(36)Φ¯j=∑i1=1N1G1i1XQji1,1∑i2=1N2G2i2XQji2,1.


## 4. 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 j, through a series of modal transformations as reported in [30], we have(37a)-λZj″+λφj′+∑J=1J0ksJ∗ZδX-XJ=wj2Zj,(37b)-1λφj″+λLr2φj-λLr2Zj′+∑J=1J0ktJ∗φδX-XJ=wj2φj.Under certain loads, we writeΦjX,τ=ZjXφjXqjτ, in modal superposition manner; (37a) and (37b) satisfy(38a)∑j=1∞Zjqj¨+wj2qj=Fp∗,(38b)∑j=1∞φjqj¨+wj2qj=Fm∗.Based on modal orthogonality, we obtain(39a)qj¨+wj2qj=Fj∗,(39b)Fj∗τ=∫01ZjFp∗dx+∫01φjFm∗dx∫01ZjX2+φjX2dx.To examine the quantitative performance, we further define the total absorbed energy of the whole system as(40)Θ=∫0T0∑j=1∞Fj∗qj˙dτ.Hence, to achieve the minimum energy absorption, (40) must yield Θ=Θmin. In virtue of the derivations in Section 3, we get Fj∗=Fj∗(τ,X1,X2,…,XN-1) and wj=wj(X1,X2,…,XN-1), all of which are the functions of (X1,X2,…,XN-1). Herein, we suppose the beam without initial disturbance and, after Duhamel integration, get the modal response as(41)qj=1wj∫0τFj∗τ=ξsin⁡wjτ-ξdξ.After substituting it to (40), we obtain(42)ΘX1,X2,…,XN-1=∑j=1∞∫0T0∫0τFj∗τFj∗ξcos⁡wjτ-ξdξdτ.Particularly for time-invariant loads, it is relatively easy to calculate the absorbed energy when system reaches into equilibrium state (w=0), which holds(43)ΘstaticX1,X2,…,XN-1=∑j=1∞Fj∗2wj2.Further, to derive the optimal layout, (42) and (43) must satisfy (44), which comes down to a minimization issue(44)∂Θ∂Xi=0,∂2Θ∂Xi2≥0.By which the optimal X1,X2,…,XN-1 can be worked out.


## 5. Numerical Verification

To demonstrate the validity of the proposed solutions and make comparisons each other, we conduct 2 typical numerical examples illustrated in Figure2 named as case 1 and case 2, respectively, in which different BCs are investigated and inherent variables are concerned, for example, the slender ratio Lr and the BC positions Xi. 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.Figure 2
Configurations of the multibeam system in case studies.


### 5.1. 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, X1=0.01 or X2=0.02, the bad condition-numbers of M or K 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 realdet⁡B=0 which is proved with pretty effectiveness. (2) “Mirror verification” (interpreted as algorithm flexibility analysis): for a set of symmetrical BCs, for example, X1=0.3, X2=0.6 and X1=0.4, X2=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.


### 5.2. Cases Studies

As shown in Figure2, 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, N=3 is employed in the study.For A-method, theB holds(1)
case1(45)B=B0I4000B112B12200B212B22200B3I41,(2)
case2(46)B=B0I1000B111B12100B211B22100B3I11.And for G-method, it holds(1)
case1(47a)G1i1=Fi1X,(47b)G2i2=∫∏j2=12X-Xj1Fi2XdX+1,(2)
case2(48a)G1i1=∏j1=12X-Xj1Fi1X,(48b)G2i2=∫X2X-12Fi2XdX+1.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 Figure3, the frequency plots of FEM and A-method at X1=0.3 and X2=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 Lr. 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.Figure 3
Frequency curves of FEM and A-method atX1=0.3, X2=0.6 condition.
(a)
Case1
(b)
Case2Figure 4
The first 12 displacemental modal shapes of A-method and FEM.
(a)
L
r
=
1 of case 1
(b)
L
r
=
1 of case 2
(c)
L
r
=
10 of case 1
(d)
L
r
=
10 of case 2
(e)
L
r
=
100 of case 1
(f)
L
r
=
100 of case 2Although the comparisons with FEM obtain satisfactory results, the high costs of A-method in computation, where the computation underN=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 (X1=0.3, X2=0.6 and X1=0.4, X2=0.7) and various N1 and N2 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 G-method. 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 KB=0, the modal sectional angles at hinge-joint BCs (X1 and X2) 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 z and the rotation φ at fix-hinge 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].Figure 5
Frequencies generated by different G-methods with variousN1 and N2.
(a)
X
1
=
0.3, X2=0.6 of case 1
(b)
X
1
=
0.4, X2=0.7 of case 1
(c)
X
1
=
0.3, X2=0.6 of case 2
(d)
X
1
=
0.4, X2=0.7 of case 2Figure 6
Modal shapes of A-method and G-methods atX1=0.3, X2=0.6 condition.
(a)
Case1
(b)
Case1
(c)
Case2
(d)
Case2Besides, it is undeniable that the slender ratio is regarded as an essential parameter in beam system [1, 32]. Hence, a large-scale variations about Lr 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 Lr<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 Lr>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 Lr, where bad computing effect occurs within low Lr, while in large Lr the consequences of G-methods are fairly acceptable. Further, the curves inherit closing regions near Lr=15 to Lr=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 Lr. On the whole, we conclude the P-method 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 Lr>30.Figure 7
Frequency curves of A-method and G-methods with change ofLr.
(a)
X
1
=
0.3, X2=0.6 of case 1
(b)
X
1
=
0.4, X2=0.7 of case 1
(c)
X
1
=
0.3, X2=0.6 of case 2
(d)
X
1
=
0.4, X2=0.7 of case 2


### 5.3. 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 Section4 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 X1 and X2, to solve the frequency is almost impossible with two unknown parameters. On the other hand, the derivations of (46) about X1 and X2 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 Figure8, we traverse X1 and X2 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’sgriddata…,“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 X1 and X2, 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 X1 and X2 is made and the two-dimension equations are conducted to get the final optimal BC values. Actually, excessive integrals must be solved at each X1 and X2 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.Figure 8
Algorithm strategy.To conduct the numeration, we setFp=10-6N/m and Fm=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 C-method, 0.02 in L-method, and 0.02 in T-method as the inherent steps for X1 or X2, respectively. Figures 9 and 10 are the frequency spectrums and the energy maps of the A-method 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.Table 1
Optimizing BCs and errors compared with A-method.
Item
A-method
P-method
C-method
L-method
T-method
X1
0.2311
0.2741
0.1983
0.182
0.2272
X2
0.7829
0.8424
0.742
0.7382
0.7296
ErrorX1
-
4.30%
−3.28%
−4.91%
−0.39%
ErrorX2
-
5.95%
−4.09%
−4.47%
−5.33%Figure 9
Frequency spectrums computed through A-method and G-methods.
(a)
A-method, initial step 0.05
(b)
P-method, initial step 0.01
(c)
C-method, initial step 0.05
(d)
L-method, initial step 0.02
(e)
T-method, initial step 0.02Figure 10
Energy maps computed through A-method and G-methods.
(a)
A-method, initial step 0.05
(b)
P-method, initial step 0.01
(c)
C-method, initial step 0.05
(d)
L-method, initial step 0.02
(e)
T-method, initial step 0.02


## 5.1. 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, X1=0.01 or X2=0.02, the bad condition-numbers of M or K 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 realdet⁡B=0 which is proved with pretty effectiveness. (2) “Mirror verification” (interpreted as algorithm flexibility analysis): for a set of symmetrical BCs, for example, X1=0.3, X2=0.6 and X1=0.4, X2=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.


## 5.2. Cases Studies

As shown in Figure2, 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, N=3 is employed in the study.For A-method, theB holds(1)
case1(45)B=B0I4000B112B12200B212B22200B3I41,(2)
case2(46)B=B0I1000B111B12100B211B22100B3I11.And for G-method, it holds(1)
case1(47a)G1i1=Fi1X,(47b)G2i2=∫∏j2=12X-Xj1Fi2XdX+1,(2)
case2(48a)G1i1=∏j1=12X-Xj1Fi1X,(48b)G2i2=∫X2X-12Fi2XdX+1.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 Figure3, the frequency plots of FEM and A-method at X1=0.3 and X2=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 Lr. 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.Figure 3
Frequency curves of FEM and A-method atX1=0.3, X2=0.6 condition.
(a)
Case1
(b)
Case2Figure 4
The first 12 displacemental modal shapes of A-method and FEM.
(a)
L
r
=
1 of case 1
(b)
L
r
=
1 of case 2
(c)
L
r
=
10 of case 1
(d)
L
r
=
10 of case 2
(e)
L
r
=
100 of case 1
(f)
L
r
=
100 of case 2Although the comparisons with FEM obtain satisfactory results, the high costs of A-method in computation, where the computation underN=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 (X1=0.3, X2=0.6 and X1=0.4, X2=0.7) and various N1 and N2 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 G-method. 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 KB=0, the modal sectional angles at hinge-joint BCs (X1 and X2) 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 z and the rotation φ at fix-hinge 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].Figure 5
Frequencies generated by different G-methods with variousN1 and N2.
(a)
X
1
=
0.3, X2=0.6 of case 1
(b)
X
1
=
0.4, X2=0.7 of case 1
(c)
X
1
=
0.3, X2=0.6 of case 2
(d)
X
1
=
0.4, X2=0.7 of case 2Figure 6
Modal shapes of A-method and G-methods atX1=0.3, X2=0.6 condition.
(a)
Case1
(b)
Case1
(c)
Case2
(d)
Case2Besides, it is undeniable that the slender ratio is regarded as an essential parameter in beam system [1, 32]. Hence, a large-scale variations about Lr 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 Lr<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 Lr>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 Lr, where bad computing effect occurs within low Lr, while in large Lr the consequences of G-methods are fairly acceptable. Further, the curves inherit closing regions near Lr=15 to Lr=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 Lr. On the whole, we conclude the P-method 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 Lr>30.Figure 7
Frequency curves of A-method and G-methods with change ofLr.
(a)
X
1
=
0.3, X2=0.6 of case 1
(b)
X
1
=
0.4, X2=0.7 of case 1
(c)
X
1
=
0.3, X2=0.6 of case 2
(d)
X
1
=
0.4, X2=0.7 of case 2


## 5.3. 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 Section4 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 X1 and X2, to solve the frequency is almost impossible with two unknown parameters. On the other hand, the derivations of (46) about X1 and X2 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 Figure8, we traverse X1 and X2 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’sgriddata…,“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 X1 and X2, 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 X1 and X2 is made and the two-dimension equations are conducted to get the final optimal BC values. Actually, excessive integrals must be solved at each X1 and X2 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.Figure 8
Algorithm strategy.To conduct the numeration, we setFp=10-6N/m and Fm=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 C-method, 0.02 in L-method, and 0.02 in T-method as the inherent steps for X1 or X2, respectively. Figures 9 and 10 are the frequency spectrums and the energy maps of the A-method 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.Table 1
Optimizing BCs and errors compared with A-method.
Item
A-method
P-method
C-method
L-method
T-method
X1
0.2311
0.2741
0.1983
0.182
0.2272
X2
0.7829
0.8424
0.742
0.7382
0.7296
ErrorX1
-
4.30%
−3.28%
−4.91%
−0.39%
ErrorX2
-
5.95%
−4.09%
−4.47%
−5.33%Figure 9
Frequency spectrums computed through A-method and G-methods.
(a)
A-method, initial step 0.05
(b)
P-method, initial step 0.01
(c)
C-method, initial step 0.05
(d)
L-method, initial step 0.02
(e)
T-method, initial step 0.02Figure 10
Energy maps computed through A-method and G-methods.
(a)
A-method, initial step 0.05
(b)
P-method, initial step 0.01
(c)
C-method, initial step 0.05
(d)
L-method, initial step 0.02
(e)
T-method, initial step 0.02


## 6. 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 ofLr>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 A-method, and the optimization is finally demonstrated to be meaningful.



---

*Source: 1269738-2018-01-16.xml*