基于Green 函数法的Timoshenko 曲梁强迫振动分析

赵 翔，周 扬，邵永波，刘 波，周 仁

(1. 西南石油大学土木工程与测绘学院，成都 610500；2. 西南油气田分公司重庆气矿开发科，重庆 400021)

我国基础设施建设的不断完善加快了土木工程、机械工程、石油工程等的发展步伐，在很多实际工程问题中，由于受场地、工程地质条件以及优化和美观设计等影响，直梁结构不再能满足工程需求，因此出现了很多曲梁结构，例如曲线桥梁、曲线隧道和弯曲机械构件等[1 − 3]。研究曲线结构的抗震与直线结构不同，它的力学特性复杂，分析起来更加困难[4]。在对这些实际工程问题的研究中，曲线结构抗震研究往往被简化为曲梁的强迫振动。许多学者已经在做这方面的研究，例如魏双科[5]建立双脊骨空间模型用于分析曲线梁桥的地震反应行为；闫磊等[6]提出了一种新型抗震体系—漂浮抗震体系，该抗震体系适用于非规则曲线桥梁的抗震；周彦良等[7]研究曲线隧道在不同地震波输入方向下不同断面的拱顶位移、内力和最大主应力变化规律；周彦良等[8]分析了曲线隧道的震害特性与机理。何燕丽、赵翔[9]建立了力电耦合的曲梁压电俘能器模型并运用Green函数求得其强迫振动的解析解。

国内外学者对曲梁的静力学研究已经比较成熟[10 − 14]，针对Timoshenko 曲梁，Lee 和Yan[15 − 16]应用位移函数法对具有强非线性边界条件的Timoshenko 曲梁的面内、外静挠度进行了研究。静力法能解决简单曲梁结构与工况的位移及内力分析问题，且可以获得解析解，但对于复杂结构、复杂工况，静力法求解难度则比较大，所以对曲梁动力特性的研究显得尤为重要。

梁的动力特性研究一直是一个经典且长期存在的问题[4, 17 − 18]。从以往的研究看，Euler、Rayleigh和Timoshenko 等各种各样的经典模型依次被提出。而曲梁的动力特性研究目前还处于探索阶段，由于曲梁结构在工程中的应用日益增多，所以越来越多的学者开始关注曲梁的振动问题，而梁的振动分为自由振动和强迫振动，其中自由振动是指振动系统按其固有频率振动，不需外力的作用，而强迫振动是指在周期性外力的持续作用下，振动系统发生的振动。目前来看针对Timoshenko曲梁强迫振动问题的研究还较少，多数的学者主要还是倾向于研究Timoshenko 曲梁的面内、外自由振动。例如，Lee[19 − 20]用伪谱法对Timoshenko曲梁进行面内、面外自由振动分析，Liu[21]基于等几何方法对Timoshenko 曲梁的面内、外自由振动进行了研究，Lv 等[22]结合改进的傅里叶级数法和瑞利-里兹法，给出了具有一般弹性边界和耦合条件的多跨Timoshenko曲梁面内振动分析的统一解。当然还有一些研究方向不同的学者，例如，Calim[23]通过对空间弯曲和扭曲的Timoshenko 梁理论公式的重写，得到了粘弹性地基上圆形梁的控制方程。

上述学者对Timoshenko 曲梁的振动问题做了很多研究，但其研究方向主要还是集中在自由振动的问题上。经过作者的文献调研，迄今为止还未有学者对Timoshenko 曲梁的强迫振动问题进行深入研究。由于结构在地震荷载作用下的破坏主要源于强迫振动，因此研究曲梁的强迫振动问题对今后的曲梁结构抗震分析是十分必要的。

从整体上看，以往的研究大多没有考虑阻尼效应。事实上，阻尼效应在工程应用中是非常重要的，Li 等[24]通过引入两个特征参数来考虑阻尼对强迫振动的影响，用Green 函数法求解了Timoshenko 直梁在强迫振动下的解析解；本文在参考文献[24]的基础上研究了具有阻尼效应的Timoshenko 曲梁(TCB)强迫振动的稳态Green 函数。从基本控制方程出发，本文依次采用分离变量法和Laplace 变换，推导出Green 函数，其中涉及的所有常数由边界条件决定。通过将某些物理量设置为零或无穷大，可以很容易地将目前的基本解简化为不存在阻尼效应的经典Timoshenko 直梁 (TB)、 Prescott 梁 (PB)和 Euler-Bernoulli 梁(EB)的基本解。通过数值计算，讨论了该方法的有效性，并给出了各种特殊几何物理量的影响。目前的解决方案以封闭和显式的形式给出，可以作为计算方法的基准。此外，根据最近的研究[25 − 27]，本文可以利用现有的解决方案对更多涉及的问题进行分析。

1 TCB 模型的建立

曲梁强迫振动的控制振动方程[28]为：

式中：v(s,t)、w(s,t)、ψ(s,t)分别为曲梁的轴向位移、径向位移和转角；p(s,t)为外部荷载；N(s,t)为轴向力；Q(s,t)为剪力；M(s,t)为弯矩；R、A、I分别表示曲梁的半径、截面面积和截面静距；E、G分别为曲梁材料的弹性模量和剪切模量；γ、µ分别表示曲梁的转动惯量和单位长度质量；κ为曲梁的剪切修正因子[29]；c1表示平动阻尼[30]；“.”表示对时间t的导数。

式中：EI和κAG分别代表弯曲刚度和剪切刚度；“'”表示对曲梁任意截面位置s的导数。假设梁上作用有如下简谐分布荷载：

将式(8a)与式(8b)联立消除Ψ，此时再将式(9)代入，即得到TCB 的振动控制方程：

2 TCB 强迫振动的稳态Green 函数

振动控制方程式(10)可以简写成以下形式：

当常数a1、a2、a3、b1、b2和b3被赋予适当的值时，Euler-Bernoulli 曲梁(ECB)和Prescott 曲梁(PCB)的振动方程也可以表示为式(11)的形式。不同梁模型所对应的a1、a2、a3、b1、b2和b3可见附录1 中的表1。

值得注意的是式(10)可以被退化为经典模型。在让c1和c2消失的情况下，设置曲梁半径R为无穷大可以得到不产生阻尼效应的传统TB 模型；如果进一步设置剪切修正因子κ为无穷大，就得到了PB 模型的控制方程；最后忽略转动惯量的影响，也就是令γ=0 可以得到EB 模型的控制方程[21]。(注*：在下文中，作者引用变量 x 代替曲梁的任意截面位置 s ，以便于与Laplace 变换中的复变量s= σ +i τ 作区分)。

根据叠加原理，式(11)的解可以表示为如下卷积积分的形式：

式中：L为梁的长度；f(x0)为外部荷载的分布函数；G(x,x0)为待求的Green 方程。从物理上来说，Green 方程G(x,x0) 指的是梁上任意一点x0作用一个单位集中力所引起的梁的响应。从数学上来说，G(x,x0)指的是下面微分方程的解：

式中，δ(·)是狄拉克函数。由式(13)可以看出G(x,x0)=W(x,x0)，即在x=x0处作用一个简谐单位集中力所引起的梁的位移W(x)。

为了推导出相应的Green 函数，作者对式(13)中的变量x施行Laplace 变换，这里引入一个新的函数We(x)：

对式(13)中的变量x施行Laplace 变换，得到：

3 确定Green 函数的待定系数

首先通过计算 ϕi(x) (i=1, 2, ···, 7)的各阶导数，来确定常数W(0)、W'(0)、W"(0)、W"'(0)、W(4)(0)和W(5)(0)，根据式(16)中Ai(x)和式(19)中 ϕi(x)的定义有：

物理上来说，式(21)建立了边界处x=0 与梁上任意截面位置x之间的关系，特别当x=L时，可以建立以下代数方程组：

如果令ϕ1(x)及其各阶导数都为0，则式(22)所描述的关系同样适用于TCB 的自由振动。

在确定了曲梁的边界条件下，常数W(0)、W'(0)、W"(0)、W"'(0)、W(4)(0)和W(5)(0)可以被求解。对不同的梁模型：EB、PB 和TCB，各种不同边界条件的表示形式参见附录2 中表2。

本文解适用于附录2 中的各种不同边界，这里以两端简支的TCB 在x=L/2 处受到简谐荷载P0cosΩt作用为例进行具体求解过程的说明，如图1所示，来说明如何确定6 个待定常数：W(0)、W'(0)、W"(0)、W"'(0)、W(4)(0)和W(5)(0)。

首先，根据x=0 处的边界条件确定位移和力学边界：

图1 两端简支的曲梁在x=L/2 处受到简谐力作用Fig. 1 A simply-supported curved beam subject to harmonic force at x=L/2

通过将位移边界式(23)代入式(8)中通过消元法得到W(0/L)、W'(0/L)、W''(0/L)、W'''(0/L)、W(4)(0/L)和W(5)(0/L)的关系式：

将其与式(21)联立可得到如下矩阵：

左边界处：x=0，此时有，W(0)=0，矩阵(25)中的前两个方程由式(24b)和式(24c)得到；右边界处：x=L，此时将式(21a)代入式(24a)，将式(21c)和式(21e)代入式(24b)，将式(21b)、式(21d)、式(21f)代入式(24c)，可依次得到矩阵式(25)中的后3 个方程。

λ1、λ2、λ3、λ4和λ5的值参见附录4，从而式(21)中的所有常数W(0)、W'(0)、W"(0)、W"'(0)、W(4)(0)和W(5)(0)被确定。因此两端简支的TCB 的Green函数的形式如下：

其中， ϕi(x)(i=1, 3, 4, 5, 6, 7)在式(19)中被给出，应该再次指出的是式(26)适用于TB、PB 和EB；通过将半径R设置为无穷大，可以简化为TB 振动模型，在此基础上，将剪切修正因子κ设置为无穷大，可以退化为PB 振动模型，最后再把转动惯量γ 设置为0，可退化为EB振动模型。

对于其它的边界条件，如两端固支，两端自由，固支-简支，悬臂等边界条件，与求解两端简支曲梁的Green 函数一样，按照同一过程也可以得到相应的Green 函数。

4 数值结果和讨论

如图1 所示,考虑到一个两端简支的TCB，梁高为h，梁长L，在x0=L/2 处受到单位简谐集中力P(x,t)=δ(x−L/2)eiΩt作用。为了便于说明，本文引入以下无量纲化参数：

4.1 Green 函数解的有效性验证

在第3 节中，已经提到TCB 模型的Green 函数解可以退化为TB 模型的Green 函数解。利用这一点可以对简支TCB模型的Green 函数解进行退化验证，看其是否与Li[24]所得到的TB 的Green函数解是一致的。因此如图2 所示，作者将本文结果的退化解与文献[24]中的简支Timoshenko 直梁模型的结果作了对照。

图2 简支TB 在中点处作用单位简谐集中力的无量纲化的挠度Fig. 2 The dimensionless deflection subject to the external unit simple harmonic concentrated force at the middle section of simply supported TB

对照显示，本文结果的退化解与文献[24]的结果完全吻合，从而本文结果的有效性得到了验证。

为了进一步验证解的有效性，利用ANSYS 软件建立了两端固支的TCB 有限元模型，并将跨中处作用单位静荷载下产生的无量纲化挠度与本研究的静态解作对比分析，如图3 所示。这里引入无量纲化参数：g1(ξ, ξ0)=G(x,x0)/wc，其中wc=L3/(192EI) 是两端固支梁的中截面x0=L/2 处受到单位集中力作用产生的最大静挠度。

图3 固支TCB 在中点处作用单位集中力的无量纲化的挠度Fig. 3 The dimensionless deflection subject to the external unit concentrated force at the middle section of fixed-fixed TCB

对比显示本文的静态解与有限元算例的位移值基本吻合，从而本文结果的有效性得到了进一步验证。

4.2 两端简支的不同梁模型共振频率研究

图4 所示是以外激励频率Ω1为自变量的无量纲化挠度g(1/2, 1/2)，图中的TB、PB、EB 均由TCB模型退化而得。

计算结果得出，当Ω1=1.0 时，EB发生共振。对于PB 和TB，它们的共振频率分别为 ΩP1B=0.997和 ΩT1B=0.984，也就是说，其一阶固有频率分别为90.138 rad/s 和88.981 rad/s，这与文献[21]中的结果一致，并且与文献[36]中所预测的值相吻合，从而进一步证明了本文所得解的正确性。

图5 所示是不同半径的TCB(L=0.5m)以外激励频率Ω1为自变量的无量纲化挠度g(1/2, 1/2)。

从图中可以看出，随着半径的增长，TCB的共振频率逐渐减小，并收敛于TB(R=∞)的共振频率。

图4 以外激励频率Ω1 为自变量的无量纲化挠度g(1/2,1/2)Fig. 4 The dimensionless deflection g(1/2, 1/2) as a function of the dimensionless frequency Ω1 of external dynamic force as independent variable

图5 不同半径的TCB 以外激励频率Ω1 为自变量的无量纲化挠度g(1/2,1/2)Fig. 5 The dimensionless deflection g(1/2, 1/2) as a function of the dimensionless frequency Ω1 of external dynamic force as independent variable for TCB with different radiuses

图6 所示是无量纲化的梁的挠曲线g(ξ, 1/2)，其横坐标是无量纲化的梁的跨度ξ=x/L。图中的TB、PB、EB 均由TCB 模型退化而得。从四种梁模型相应Green 函数曲线的模态形状可以看出，Ω1=0.5 可以激发一阶模态。在该频率激励下，PB与EB 的解极为接近，这说明转动惯量对Green 函数的影响较为微弱。相反的，TB 与EB 的解差异显著，这说明剪切效应对梁的振动挠曲线影响较大。如图6 所示，TCB 与三种直梁模型的挠曲线相似，其挠度比TB 模型的挠度小。但是从下一小节的研究结果中不难发现随着半径R的变小，TCB 的挠曲线形状会发生较大的变化。

4.3 不同半径下简支TCB 的挠度变化

图7 所示是不同半径的TCB 无量纲化的挠曲线g(ξ, 1/2)，其横坐标是无量纲化的梁的跨度ξ=x/L。

图6 在外激频Ω1=0.5 作用下，EB、PB、TB 和TCB(R=5)的无量纲化挠度g(ξ,1/2)Fig. 6 The dimensionless deflection g(ξ,1/2) of EB, PB, TB,and TCB (R=5) corresponding to Ω1=0.5 of external dynamic force.

图7 不同半径下，TCB 的无量纲化挠度g(ξ,1/2)Fig. 7 The dimensionless deflection g(ξ,1/2) of TCB with different radiuses

从图中可以看出，随着半径R的减小，TCB的挠度逐渐变小，应当特别指出的是：当半径R=0.25 m 为梁长L=0.5 m 的1/2 时，其边界处转角为0，挠曲线的形状类似于两端固支的直梁。

图8 所示是以半径R为自变量的无量纲化挠度g(1/2, 1/2)。

图中可以看出随着半径R的增长，在起始阶段TCB 的跨中挠度呈线性增加，随后增加趋势逐渐变得缓和，最终收敛于相同条件下TB 模型的无量纲化挠度。

4.4 阻尼效应对简支TCB 的影响

图9 与图10 所示是以半径R=2 m，梁长L=0.5 m的TCB 为例，分别以阻尼比ζ1和ζ2为自变量的无量纲化位移g(1/2, 1/2)。

图8 以半径R 为自变量的TCB 的无量纲化挠度g(1/2,1/2)Fig. 8 The dimensionless deflection g(1/2,1/2) of TCB with radius R as independent variable

图10 以阻尼比ζ2 为自变量的无量纲化挠度g(1/2, 1/2)Fig. 10 The dimensionless deflection g(1/2,1/2) of TCB with damping ratio ζ2 as independent variable

从图中可以看出：随着阻尼比ζ1和ζ2的增加，位移值变得越来越小，这是符合客观物理实际的；另外，通过对比图9 与图10 还可以看出平动阻尼的影响比转动阻尼要大。控制方程中的平动阻尼代表空气阻尼，而转动阻尼代表材料阻尼[33]。计算结果表明，材料阻尼对简支曲梁的位移响应影响较小，因此可以忽略。

5 结论

本文通过Laplace 变换，系统地研究了两端简支、固支、自由边界条件下TCB 强迫振动的Green函数。并以两端简支边界条件为例，证实了所得到的TCB 强迫振动Green 函数可以退化到TB、EB和PB 的强迫振动Green 函数。并通过数值手段，研究了以上四种梁模型的共振频率、不同半径下TCB 的挠度变化和阻尼效应的影响。并得出以下几点重要结论：

(2)随着半径R的增大，TCB 的无量纲挠度逐渐变大，并最终收敛于TB 的无量纲挠度，应当特别指出的是：当半径R=0.25 m 为梁长L=0.5 m 的1/2 时，其边界处转角为0，挠曲线的形状类似于两端固支的直梁。

(3)随着阻尼比ζ1和ζ2的增加，TCB 的跨中挠度近似于线性越小。此外，由于材料阻尼对简支曲梁的位移响应影响较小，因此可以忽略。

在解的验证部分，所得到的TCB 的Green 函数退化为TB 的Green 函数后与Li 等[24]所得到的TB 的Green 函数完全吻合，且退化为EB 的Green函数与Abu-Hilal[33]所得到的EB 的Green 函数一致。此外，本文的静态解与有限元算例的位移值基本吻合，从而本文结果的有效性得到了进一步验证。

本文的研究结果可为相关理论研究和工程应用提供有益参考，文中所涉及的物理量符号说明参见附录3 中表3。