相比较于拉普拉斯变换，对于复杂型式的载荷，例如移动载荷，其数学表述与空间变量存在耦合特性，使得拉普拉斯反变换的解析解很难求得，即使在固定的阶跃载荷情况下也得运用复杂的谱分析法才能获得解析解[9]。且对于数值拉普拉斯反变换，例如Dubner和Durbin法[14-15]，对于参数a和T的选择也较依赖经验公式，且效果不一定好。
1 建模与求解 1.1 建模 对于二自由度铁木辛柯梁，截面弯曲引起的位移zb和剪切引起的位移zs可以表示为
 | (1) |
式中：E和G分别为梁的弹性模量和剪切模量；I为截面惯性矩；ρ为材料密度；A为截面积；k为剪切系数[16], 对于矩形截面，一般取5/6；zb和zs共同构成梁的总位移ZT；x为梁的横向坐标；t为时间。为便于推导，引入截面转动角φ和剪切角γ。
进一步，考虑K-V阻尼模型[17]:
 | (2) |
式中：σ和ε分别为梁内部的拉伸应力和应变；τs为梁内部的剪切应力；Cb和Cs分别为梁的拉伸和剪切黏性阻尼系数。采用截面线性应变模型：
 | (3) |
式中：h为截面上任意一点到中线的距离，为积分变量。借此，截面弯矩Mb和剪力Q可以显式地表示为
 | (4) |
式中：B为梁截面y向总长；H为z向总长。
因此，考虑K-V阻尼模型的铁木辛柯梁动力学方程可以表示为[13]
 | (5) |
式中：F为外力；q=[zs, zb]T。
式(5)同时考虑了外力F，可以包括任意型式的力和力矩，但本文仅考虑阶跃载荷和移动载荷。因此，F=[f(x, t), 0]T。为方便研究，将式(5)进行无量纲化，引入无量纲化的各变量：
因此，式(5)可以改写为
 | (6) |
式中：F*为F的无量纲形式。
1.2 求解 定义
 | (7) |
先验证系统是否具有比例阻尼特性，分别对KM-1C和CM-1K进行求解，分别为
 | (8) |
 | (9) |
可以得出，当Cs*/Cb*=λ02成立时，系统具有比例阻尼特性，即，KM-1C=CM-1K。因此，为了研究阻尼对系统的总体影响，这一假设可以大大降低求解的难度，而且也可以得到实用的结果。Chen在分析系统参数对轴向受扭状态下的铁木辛柯梁的频率特性研究中，将Cb和Cs处理为相等的常数[11]。而Capsoni等证明具有临界振动特性的Cb和Cs的大小与系统的参数密切相关，它们之间比值关系会随着参数发生变化[18]。因此，本文的假设不会使系统的机理模型产生矛盾。
1.2.1 无阻尼振动求解 本节之后的推导仍然同时保留Cb和Cs，但其已满足比例阻尼条件，假设式(7)无阻尼的解为[Zs, Zb]T=[C1, C2]Teiw0ζ。因此，根据非平凡解[C1, C2]的存在性，可得特征方程：
 | (10) |
式中：Ω=w2；λ为特征方程式变量，可得
 | (11) |
式(10)解的结构决定其具体形式。一般的，其通解为
 | (12) |
式中：ai为系数；Zi的下标i为s或b；λ1和λ2分别为式(11)的2个正根。
将其转化至实域时，可由sinh(?)、cosh(?)、sin(?)和cos(?)等实函数进行表示。本文以两边固支的边界条件为例，其他形式的边界条件可以用同样的方法得出。其中，固支边界条件为
 | (13) |
式中：ZL为z向总位移。
1) 当Ω＜λ0Lr2时，式(10)由一对实根和一对共轭复根构成，分别记为±η和±α。其频率方程和特征向量分别为
 | (14) |
 | (15) |
式中：
2) 当Ω=λ0Lr2时，式(10)由2个0重根和一对虚根±α组成。其频率方程和特征向量分别为
 | (16) |
 | (17) |
式中：
3) 如果Ω＞λ0Lr2，式(10)的解由一对虚根±α和另一对虚根±β构成，其频率方程和特征向量分别为
 | (18) |
 | (19) |
式中：
需要指出的是，考虑到实际情形，式(17)的情况不存在[6, 8]。
1.2.2 阻尼振动求解 同1.2.1节，根据分离变量法，比例阻尼条件下的系统模态响应为
 | (20) |
式中：Qn(ζ)为比例阻尼情形下的模态响应函数，仅和时间变量ζ相关。同时，考虑到式(6)中的外力源项，此处将式(20)代入式(7)，然后在等式两侧同时左乘[Zsm, Zbm](m≠n，m、n为正整数)，并将所得的结果对X从0到1积分，根据模态振型的正交特性，式(7)最终可以被转换为第m阶模态的振动表达式：
 | (21) |
其中：
 |
此处，式(21)的显式结果可以通过Zhang等的求解过程导出[13]。其中，利用Duhamel积分或者卷积积分法，系统在受载时的响应可以表示为如下结果：
1) 阻尼比满足κm=0时,
 | (22) |
2) 阻尼比满足0＜κm＜1时,
 | (23) |
式中：上标d代表阻尼情形。
3) 阻尼比满足κm=1时,
 | (24) |
4) 阻尼比满足κm＞1时,
 | (25) |
具体的，位于X=d*阶跃载荷可以通过海威塞得函数H(ζ)定义为
 | (26) |
移动载荷可以定义为
 | (27) |
2 算例分析 图 1是关于Durbin法的一个数值实验，可以看出当乘积aT一定时，随着T的改变，梁内部物理量求解误差相差极大。而模态叠加法，因其具有较好的数值稳定性，常用作检验数值方法的标准和较强的适用性，因此被第1节采用。
 |
| 图 1 Durbin数值反拉普拉斯求解效果示意图 Fig. 1 Schematic of solving accuracy of Durbin's inverse Laplace |
| 图选项 |
2.1 模态求解 本节首先对频率和模态特性进行计算，参考一组真实原始数据[13]：L=1 800 m，A=20 m2，I=309.156 3 kg·m2，E=3.63×1010 N·m-2，G=1.578 3×1010 N·m-2，ρ=5.128 2×103 kg·m-3，Cb=2×108 Pa，k=5/6，F0=2×104 N。
图 2(a)、(c)为在无阻尼时的频率求解结果，随着Lr的增大，系统的固有频率逐渐降低；此外，从图 2(b)、(d)的阻尼变化趋势可以看出：采用铁木辛柯梁模型计算出的阻尼比，在Lr=10时，系统在前3阶频率内为振荡阻尼特性，在其之后进入了过阻尼状态，并随着阶次的增长迅速增大；而当Lr=100时，系统在第48阶模态之后才进入过阻尼振动状态；而且，与Lr=10同阶模态相比，Lr=100时的阻尼值明显减小。这说明系统在受载后将更慢恢复到平衡状态。对于所研究的2种长细比，过阻尼状态均发生在第1种和第3种情况的交界处，当具有第1种情况的振型特性时，过阻尼状态不会发生。
 |
| 图 2 模态频率和阻尼比变化趋势图 Fig. 2 Changing tendency of modal frequencies and damping ratios |
| 图选项 |
图 3(a)~(f)给出了系统的模态振型图。其中，图 3(a)~(c)和图 3(d)~(f)分别为在Lr=10和Lr=100时的模态振型，图 3(a)、(b)、(d)、(e)为Ω＜λ0Lr2时的解，图 3(c)、(f)对应Ω＞λ0Lr2的情况。从模态阵型形状来看，第1类振型呈现出正余弦波动的特性，第3类振型则呈现出2种不同正弦叠加的特性。对于第1种情况，阵型函数更为平滑。
 |
| 图 3 归一化模态振型 Fig. 3 Normalized modal shapes |
| 图选项 |
2.2 阶跃载荷 本节对系统在承受阶跃载荷情况下的动态响应进行了研究，并对Lr=10和Lr=100的2种不同长细比进行对比分析。选用X=0.2和X=0.8为2个不同的研究点，分别对对称受载d*=0.5和非对称载荷d*=0.25进行了探究。其中，系统的总响应为前500阶模态叠加的结果。
从图 4、图 5中可以看出，在阶跃载荷的作用下，阶次越高，模态响应所占的比例越少，证明了在非简谐特性的阶跃载荷作用下，系统表现出了低频模态为主导的特性。从模态响应的变化趋势来看，由于系统在高Lr时具有较低的阻尼比，其在受载后回至平衡状态的时间也大大增加。Lr=10和Lr=100的模态响应趋势也验证了这一特点。系统在对称阶跃载荷的作用下，其响应相对于中点，理论上也应为对称的，如图 6所示，X=0.2和X=0.8的求解结果证实了这一结论。而当载荷被偏置后(d*=d/L=0.25)，从图 7可以看出系统的响应随时间的变化发生了变化，然而2个研究点处的振幅具有相同的量级。
 |
| 图 4 Lr=10模态响应 Fig. 4 Modal response when Lr=10 |
| 图选项 |
 |
| 图 5 Lr=100模态响应 Fig. 5 Modal response when Lr=100 |
| 图选项 |
 |
| 图 6 d*=0.5时无量纲位移总响应 Fig. 6 Non-dimensionalized response of total displacement when d*=0.5 |
| 图选项 |
 |
| 图 7 d*=0.25时无量纲位移总响应 Fig. 7 Non-dimensionalized response of total displacement when d*=0.25 |
| 图选项 |
此外，通过模态响应图或总响应图，对比无阻尼振动和有阻尼振动的结果可以得出：当Lr较小时，系统的各阶阻尼系数较大，振动收敛很快，例如对于Lr=10，响应在3~4个波峰之后基本趋于稳定，而对于大Lr情形，系统则需要历经更长的时间才能收敛。这也间接证明了对于阶跃响应，系统的低阶模态起主导作用。
2.3 移动载荷 本节对系统承受移动载荷时进行了进一步的研究。载荷的计算条件为v0=0.001 m/s。计算时间考虑载荷沿起始位置X=0运动至终点X=1。
图 8和图 9的计算结果表明，无论在低Lr或高Lr情况下，无阻尼响应曲线在移动受载时均呈现出了往复波动的特性。而得益于在低Lr时系统具有高阻尼比，阻尼响应在Lr=10时没有呈现出波动特性。对于高Lr，从响应结果可以看出，较低的模态阻尼比使得系统在考虑阻尼的情况下也呈现出明显的波动特性。而且，当载荷运动至中间位置和末端位置时，其响应竟远远超过了无阻尼振动的响应，这是因为当载荷在低速运动时，阻尼的存在使得考虑阻尼振动的梁频率比自然振动频率更低，载荷的低频分量更接近有阻尼振动时的频率，从而导致局部更接近共振频率所致。
 |
| 图 8 Lr=10无量纲位移总响应 Fig. 8 Non-dimensionalized response of total displacement when Lr=10 |
| 图选项 |
 |
| 图 9 Lr=100无量纲位移总响应 Fig. 9 Non-dimensionalized response of total displacement when Lr=100 |
| 图选项 |
3 结论 本文采用了实模态叠加法导出了铁木辛柯梁在承受阶跃载荷和移动载荷时的解析解，其中，通过比例阻尼条件将K-V模型转化为了比例阻尼模型。对动态响应的分析表明：
1) 铁木辛柯梁模型可以预测更复杂的模态振型形状，因而其频率求解更准确。
2) 系统在承受阶跃载荷时，系统的响应主要由低阶频率主导，对于小长细比，阻尼振动可以很快收敛，而对于大长细比，振动收敛很慢，这是因为当长细比较大时，各阶模态的阻尼比明显减小。
3) 系统在承受低速移动载荷时，无阻尼振动和有阻尼振动的结果表明，小长细比情况下，有阻尼振动由于能量很快被耗散掉，所以响应没有呈现出波动特性，而在大长细比下，有阻尼振动波动十分明显，而且小速度移动的载荷引起了有阻尼振动的共振，使得振幅远远超过了无阻尼的响应。
参考文献
| [1] | LABUSCHAGNE A, RENSBURG N F J V, MERWE A J V D. Comparison of linear beam theories[J].Mathematical & Computer Modelling, 2009, 49(1-2): 20–30. |
| [2] | STEPHEN N G. The second spectrum of Timoshenko beam theory-Further assessment[J].Journal of Sound and Vibration, 2006, 292(1-2): 372–389.DOI:10.1016/j.jsv.2005.08.003 |
| [3] | HAN S M, BENAROYA H, WEI T. Dynamics of transversely vibrating beams using four engineering theories[J].Journal of Sound and Vibration, 1999, 225(5): 935–988.DOI:10.1006/jsvi.1999.2257 |
| [4] | KOCATVRK T, ?IM ?EK M. Dynamic analysis of eccentrically prestressed viscoelastic Timoshenko beams under a moving harmonic load[J].Computers & Structures, 2006, 84(31): 2113–2127. |
| [5] | BOLEY B, CHAO C. Some solutions of the timoshenko beam equations[J].Journal of Applied Mechanics, 1955, 25: 579–586. |
| [6] | XIN Y F, SONG Y C, ZHU D C, et al. Elastic impact on finite Timoshenko beam[J].Acta Mechanica Sinica (English Series), 2002, 18(3): 252–263.DOI:10.1007/BF02487953 |
| [7] | 邢誉峰. 有限长Timoshenko梁弹性碰撞接触瞬间的动态特性[J].力学学报, 1999, 31(1): 68–74. XIN Y F. The characteristics of Timoshenko beam during the process of elastic impact and contact[J].Acta Mechanica Sinica, 1999, 31(1): 68–74.(in Chinese) |
| [8] | SU Y C, MA C C. Transient wave analysis of a cantilever Timoshenko beam subjected to impact loading by Laplace transform and normal mode methods[J].International Journal of Solids & Structures, 2012, 49(9): 1158–1176. |
| [9] | SU Y C, MA C C. Theoretical analysis of transient waves in a simply-supported Timoshenko beam by ray and normal mode methods[J].International Journal of Solids & Structures, 2011, 48(3): 535–552. |
| [10] | HU M Y, WANG A W, ZHANG X M. Approximate analytical solutions and experimental analysis for transient response of constrained damping cantilever beam[J].Applied Mathematics and Mechanics, 2010, 31(11): 1359–1370.DOI:10.1007/s10483-010-1368-9 |
| [11] | CHEN W R. Parametric studies on bending vibration of axially-loaded twisted Timoshenko beams with locally distributed Kelvin-Voigt damping[J].International Journal of Mechanical Sciences, 2014, 88: 61–70.DOI:10.1016/j.ijmecsci.2014.07.006 |
| [12] | GU L L, QIN Z Y, CHU F L. Analytical analysis of the thermal effect on vibrations of a damped Timoshenko beam[J].Mechanical Systems & Signal Processing, 2014, 60-61: 619–643. |
| [13] | ZHANG X Y, ZHU M, LIANG H Q. Dynamic analysis of the continuous fluid-structure system based on Timoshenko model and considering damping: AIAA-2017-1984[R]. Reston: AIAA, 2017.https://www.researchgate.net/publication/313455622_Dynamic_analysis_of_the_continuous_fluid-structure_system_based_on_Timoshenko_model_and_considering_damping |
| [14] | DUBNER H, ABATE J. Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform[J].Journal of the ACM, 1968, 15(1): 115–123.DOI:10.1145/321439.321446 |
| [15] | DURBIN F. Numerical inversion of Laplace transforms:An efficient improvement to Dubner and Abate's method[J].Computer Journal, 1974, 17(4): 371–376.DOI:10.1093/comjnl/17.4.371 |
| [16] | JENSEN J J. On the shear coefficient in Timoshenko's beam theory[J].Journal of Applied Mechanics, 1966, 33(2): 621–635. |
| [17] | ZHAO H L, LIU K S, ZHANG C G. Stability for the Timoshenko beam system with local Kelvin-Voigt damping[J].Acta Mathematica Sinica(English Series), 2005, 21(3): 655–666.DOI:10.1007/s10114-003-0256-4 |
| [18] | CAPSONI A, VIGANòG M, BANI-HANI K. On damping effects in Timoshenko beams[J].International Journal of Mechanical Sciences, 2013, 73(8): 27–39. |
