受单点横向非定常约束梁的响应分析

丁维高，谢 进

(西南交通大学 机械工程学院，成都 610031)

梁在不同边界条件约束下的响应问题是振动理论研究领域的经典问题。在力学中，将显含时间t的几何约束称为非定常约束，否则为定常约束[1]。受定常约束梁在各种激励下的响应分析已经有了大量研究成果。简单的边界条件如铰支、固支、自由[2-5]。复杂的边界条件如多跨梁[6-7]、弹性支撑[8-10]等。

但在工程中存在许多梁受非定常约束的工况。例如梁式俘能器被振动体所驱动时，可以认为梁的夹紧端受到了振动体的非定常约束[11]；在机械运动学的研究中，原动件如果是由控制电机驱动的，则机构的运动也受非定常约束，比较典型的机械装置是由伺服系统直接驱动和控制的铁路轨道梁转辙机，轨道梁便受到了转辙机的非定常约束作用[12-13]；在柔顺机构中，也存在柔顺杆上某特定点的运动为已知时间函数的工况[14]。

尽管在时变结构动力学[15-16]中，对梁的伸展问题等非定常边界问题已有一些研究成果[17-18]，但对于受横向非定常约束梁的研究还较少，且多为非定常约束作用于梁的端点时的工况。在这种情况下，通常的处理方法是将问题视为为经典“承受支撑运动的梁”问题[19]，再使用位移影响函数法[20]求解。然而，对于非定常约束作用在梁上其他任意位置情况，尽管可以先将模型转化为多跨梁问题，再利用位移影响函数的求解，但求解过程较为繁琐。因此，本文提出：直接使用第一类拉格朗日方程来解决该问题，进而使用简单边界条件的模态函数来对问题进行求解和分析。

本文首先利用第一类拉格朗日方程与欧拉-伯努利梁理论建立梁动力学方程，在此情况下，系统的动力学方程为微分代数方程(Differential-Algebraic Equation, DAE)[21-22]。如果不考虑非线性因素，单点非定常约束下梁的代数方程与微分方程均为线性方程。本文将利用求解线性强迫振动常微分方程(Ordinary Differential Equation, ODE)稳态解的求解方法求解出这类线性DAE问题的解析解：①考虑到工程中大量的结构振动或机械运动具有周期性，而周期性函数可通过三角级数将其展开为多个谐波函数来表示，故本文研究梁所受到的非定常约束为谐波函数的情况；②对解析解与数值解进行比较，并作为特例应用于求解承受支撑运动的梁，以验证解析解的正确性及普适性；③对梁非定常约束作用点对于梁的动力学性能的影响进行分析和讨论。

1 受单点横向非定常约束梁的动力学模型

1.1 单点横向非定常约束下梁动力学方程的建立

图1为受单点非定常位移约束梁的计算模型。长度为l的梁端点处自由或被定常约束限制，同时，在梁的长度方向s=al，0≤a≤1处受有非定常约束。非定常约束表示为s处的横向位移为随时间变化的函数P(t)。假设P(t)的幅值远小于梁的长度，则可以忽略梁的非线性因素的影响，梁的横向运动可以用模态函数与其对应的模态坐标表示。

图1 受单点非定常约束梁的计算模型Fig.1 Calculation model of beam under one-point transverse rheonomicis

设仅考虑定常约束时梁的第i阶模态函数为φi(x)，第i个模态坐标响应为qi(t)，则非定常约束方程可表示为

(1)

均质欧拉-伯努利梁模型的振动偏微分方程为

(2)

式中：m=ρA为单位长度梁的质量，ρ为密度，A为截面面积；E为弹性模量；I为极惯性矩；cs为线性材料阻尼系数；ca为空气阻尼系数。式(2)的通解可以表示为

(3)

第i阶模态函数归一化后的一般形式为

(4)

式中：l为梁长度，λi为梁特征方程的解，C1,i～C4,i仅与定常边界条件和特征值λi有关。梁的弹性势能可表示为

(5)

梁的动能可以表示为

(6)

阻尼力的虚功可以表示为

(7)

根据第一类拉格朗日方程

(8)

可得到受单点非定常约束梁模态坐标下的动力学方程

(9)

式中：κ(t)为拉格朗日乘子；ωr为无阻尼固有频率

(10)

ζr为模态阻尼比

(11)

1.2 谐波约束下动力学方程的解析解

由于工程中大量的结构振动或机械运动均具有周期性，故本文研究中假设约束函数P(t)为一个周期性连续有界函数：P(t)=flcos(ωt)，其中f为一个无量纲系数。

在此假设下，可以期望形如式(9)的线性DAE方程具有周期性的稳态解析解，并且稳态解也为谐波函数。

仿照求解强迫振动ODE稳态响应的方法[23],假设第r阶模态坐标稳态解为

qr(t)=Arsin(ωt)+Brcos(ωt)

(12)

由于非定常约束是谐波函数，则拉格朗日乘子也具有谐波形式

κ(t)=κAsin(ωt)+κBcos(ωt)

(13)

将式(12)与式(13)代入式(9)中得到

(14)

(15)

比较式(14)与式(15)中两边的sin(ωt)和cos(ωt)的对应系数，得到系数的线性方程组

(16)

(17)

(18)

(19)

由式(16)、式(17)可以使用分离出模态坐标系数

(20)

(21)

将式(20)、式(21)代入式(18)、式(19)中得到

(22)

(23)

式(22)、式(23)是关于κA，κB的线性方程，解得

(24)

(25)

将式(24)与式(25)代回式(20)、式(21)即可得到Ar与Br。假设激励频率ω与第一阶、不为0的频率激励频率之比为1/β

(26)

梁的稳态响应函数可以表示为

ω313,r(β,ζr)cos(ωt)]

(27)

式(27)中ψr(al)，ω212,r(β,ζr)，ω313,r(β,ζr)具体表达式见附录A。需要特别指出的是：式中质量m被约去，梁的动态响应的仅与约束位移幅值fl、模态阻尼ζr、一阶模态频率比1/β有关，且与位移约束幅值呈线性正比关系。这与受简谐力激励的梁的性质类似。

式(11)所确定的阻尼系数及各阶模态阻尼可表示为一阶、二阶模态阻尼ζ1，ζ2的函数

(28)

特别的，当忽略黏性空气阻尼时，各阶模态阻尼可表示为一阶模态阻尼ζ1的函数

(29)

总之，如果确定了约束位置幅值fl与作用位置al，可以通过对一阶频率比1/β和一阶、二阶模态阻尼ζ1，ζ2的变化来研究梁的模态响应的变化，特别是当忽略空气阻力时，仅需讨论一阶频率比1/β和一阶模态阻尼ζ1。

为了便于研究，在后续的计算中使用式(29)来计算梁的各阶阻尼的变化。

1.3 解析解数值验证及其在支撑运动梁中的应用

通过将1.2节中得到的解析解与数值解的对比，以及将经典的承受支撑运动的梁作为本文研究的受有单点非定常约束的解析解的应用特例，说明本文所得到的解析解的正确性及普适性。

1.3.1 解析解的数值验证

以一端固定的悬臂梁为例，使用4阶龙格库塔法求解式(9)，指标约简是将代数方程对时间求导的实现的。取梁的参数l=1 m，m=1 kg/m，f=0.1，a=0.5，ζ1=0.1,β=2，ω=1 rad/s。数值计算时取模态的数量为前7阶。

悬臂梁的模态函数表达式为

(30)

将参数与式(30)代入解析解，式(27)中得到

w(x,t)=[-1.353 0cos(t)-0.003 5sin(t)]ψ1(x)+

[0.054 8cos(t)-0.001 6sin(t)]ψ2(x)+

[0.000 2cos(t)]ψ3(x)

……

(31)

其中特征值由

cosλcoshλ+1=0

(32)

求出。图2(a)、图2(b)分别给出了一阶、二阶模态响应数值计算结果与解析解的对比。一阶模态响应的解析解与数值解之间的差值变化如图3所示。

图2 模态时域响应的解析解与数值解比较Fig.2 Contrast diagram of time domain response conclude by analytic solution and numerical solution

由于数值解包括了梁的瞬态响应部分，则开始部分解析解与数值解之间的差别较大，经过5 s之后两者误差的数量级会小于10-7。由于在DAE数值求解中存在有约束违约现象[24]，所以，经过较长时间的迭代后，两者差的绝对值会有增大的趋势。然而，这些并不妨碍得到解析解与数值解基本一致的结论，也并不妨碍利用解析解进行梁的动态响应分析。

图3 一阶模态时域响应解析解与数值解的误差 Fig.3 Error of time domain response of mode 1 conclude by analytic solution and numerical solution

1.3.2 在支撑运动梁中的应用

经典的“承受支撑运动的梁”问题，如图4所示。对于此问题，可采用位移影响函数法求解。Erturk等表明：位移影响函数法得到的解与实验吻合，可以准确地预测梁的动态响应。图4的计算模型可以视为本文研究的受有单点非定常约束、当约束作用点在梁的端点时的特例。利用本文提出的方法，可以得到图4的梁解析解。梁的定常边界条件为左端转角固定，右端自由

(33)

图4 承受支撑运动的梁的计算模型Fig.4 Calculation model of beam under support campaign

由x=0处的边界条件，得到C4,i=C3,i=0，则特征方程为

cosh(λ)sin(λ)+cos(λ)sinh(λ)=0

(34)

数值求解式(34)得到λ=0，2.365 0，5.497 8，8.639 4，11.781 0，14.922 6，18.064 2，…。

归一化的第i阶模态函数为

(35)

特别是，对于0特征值λ0=0，模态函数为

(36)

非定常约束为

(37)

由悬臂梁特征方程，式(32)可以得到悬臂梁的特征值为λ*=1.875 1，4.694 1，…。

图5 梁末端点位移无量纲频率响应Fig.5 Dimensionless displacement response of the end point on beam

2 谐波非定常约束与激励力作用的差异

本节对谐波非定常约束与谐波外激励对梁的作用效果进行比较。悬臂梁的右端被一个刚性正弦机构所带动，曲柄O5O6匀速转动，其长度足够短，以至于不需要考虑梁的几何非线性效应。点O2的位移为曲柄O5O6转角的正弦函数，为简便起见，称此情况为“约束工况”，如图6(a)所示；悬臂梁在右端受到一个谐波周期力F(t)=pEIcos(ωt)激励的作用，其中p为无量纲常数，称此情况为“力工况”，如图6(b)所示。

图6 悬臂梁受谐波非定常约束Fig.6 Cantilever beam under transverse rheonomicis

根据经典的振动理论，力工况的第r阶模态坐标方程为

(38)

其稳态响应解为

(39)

首先考虑a=1的简单情况，此时力和非定常约束都作用于悬臂梁的末端。计算过程中模态数量取为7,取f=0.1，取p=0.001，ω=1。考虑工程上最为常见的欠阻尼情形，取ζ1=0.001，则ζ2=0.063，ζ2=0.175，此时前三阶模态阻尼均为欠阻尼。图7给出了两种工况下前三阶响应的对比。

图7 ζ1=0.001，a=1时约束工况和力工况的前三阶模态响应Fig.7 First three orders model response under transverse rheonomicis and force(ζ1=0.001，a=1)

当频率比为1时，力工况的一阶模态响应达到了峰值，而约束工况的一阶响应并没有一个明显的峰值，二阶、三阶模态响应却达到了一个极小值；力工况的二阶模态达到峰值时，约束工况的一阶和三阶模态达到极小值；力工况的三阶模态达到峰值时，约束工况的一阶和二阶模态达到极小值。而在约束工况中，一阶、二阶和三阶模态同时达到峰值，并且达到峰值时的频率比处于力工况的两阶共振频率比之间。非定常约束对模态响应的作用可以解释为：若激励频率达到某阶共振频率，对应此阶数的模态响应最为明显，而由于弹性体的运动受到了非定常约束的限制，其余阶数的模态则必须降低响应幅值，以提升其在共振模态的幅值中所占的比例。进而可知，在激励达到某阶模态的主频时，其余阶模态会产生一个极小值。如果某阶模态响应达到了峰值，其余阶模态为满足约束，也必须进行相应的补偿。

3 约束作用位置对悬臂梁动态响应的影响

3.1 共振频率位置的估算及解曲线

当非定常约束在梁的任意位置时，如何计算共振频率的位置，以及用解曲线表示梁的动态响应。

图8所示为ζ1=0.001，a=1时梁上0.3l，0.5l和0.7l处三个位置的频率响应幅值曲线。将图8与图7中的约束工况进行比较可知：当各阶模态达到共同的峰值时，梁上各点的响应幅值也会达到响应的峰值。因此，可以认为在模态的共同峰值处发生了共振。

假设第r阶模态频率响应幅值为Hr。由式(27)有：

(40)

[-ω32(β,ζr)ω13,r(β,ζr)+ω31(β,ζr)ω23,r(β,ζr)]2=

[ω31(β,ζr)2+ω32(β,ζr)2][ω13,r(β,ζr)2+ω23,r(β,ζr)2]

(41)

设

(42)

将式(42)代入式(41)中得

(43)

γ′4(β)=γ2(β)[2γ′2(β)γ1(β)+γ′1(β)]+
γ3(β)[2γ′3(β)γ1(β)+γ′1(β)]=0

(44)

取a=1，ζ1=0.001，通过二分法求式(44)的数值解，得到1/β={1.000,4.388,6.267,14.224,…}，分别对应图7约束工况中模态响应的从左到右的第一个极小值，第一个峰值，第二个极小值，…。

图8 ζ1=0.001，a=1时梁上0.3l、0.5l、0.7l处的位移响应幅值Fig.8 Displacement response of beam on 0.3l，0.5l，0.7l (ζ1=0.001，a=1)

注意到当各阶模态都处于欠阻尼情况下时，ζ1的值会非常小，可以认为γ3(β)很小，而ψr(al)是有界的。因此可以近似的认为：为使式(44)成立，仅需

γ2(β)=0

(45)

取a=1，ζ1=0.001，由式(45)得到1/β={1.000,4.388,6.266,14.254,…}，与二分法求得的数值解非常相近。由此可知，在欠阻尼情况下，可以使用式(45)计算梁响应的峰值频率比。

为了研究非定常约束作用在不同位置时梁的响应，首先选定一阶模态阻尼，然后遍取非定常约束的作用位置参数a在区间(0,1)中的值。取定a的值，设定1/β的搜索范围以求出满足式(45)的1/β值。以a值的作为横坐标，搜索得到的1/β的值作为纵坐标，就可以得到表示不同约束作用位置下响应为峰值和极小值时的频率比的值的曲线，简称为“解曲线”。ζ1=0.001，1/β∈(0,20)的解曲线如图9所示。图9中粗实线为响应极小值，细实线为响应峰值。

图9 ζ1=0.001，1/β∈(0,20)的解曲线Fig.9 Graphs of solution(ζ1=0.001，1/β∈(0,20))

从图9可知，各个响应极小值频率比随参数a的不同变化很小，并且总是在梁的固有频率附近；而响应峰值频率比受参数a影响变化很大，但各个响应峰值总在梁的两阶相邻固有频率附近变化。响应峰值的分布规律符合本征频率的Rayleigh-Courant-Fisher定理[25]

另外，解曲线中有一些间断的位置，间断的位置在使ψr(al)=0的a*附近。例如ψ2(al)=0的解为a*=0.783 4，…。为分析间断所代表的物理意义，在间断a*位置的两侧取a={0.7,a*,0.8, 0.9}分别计算出前三阶模态的频率响应，如图10所示。当a=0.7时，在图9中穿过了5次解曲线，对应着频率比从0增加到20的过程中，图10(a)中的三个极小值和两个峰值的位置;当a=0.8和a=0.9时，在图9中一次穿过间断与三次穿过解曲线,穿过间断的位置分别对应于图10(c)中1/β=6.3附近与图10(d)中1/β=17.5附近。此两处近乎重合的峰值说明了在此位置时峰值和极小值趋于合并。

当a=0.7时，在图9中穿过了5次解曲线，对应着频率比从0增加到20的过程中，图10(a)中的三个极小值和两个峰值的位置;当a=0.8和a=0.9时，在图9中一次穿过间断与三次穿过解曲线,穿过间断的位置分别对应于图10(c)中1/β=6.3附近与图10(d)中1/β=17.5附近。此两处近乎重合的峰值说明了在此位置时峰值和极小值趋于合并。实际上，ψr(al)=0的第r阶模态会在即将合并的频率处有一个较为明显的峰值，当a=a*时，如图10(b)所示。峰值与极小值完全合并在了一起。此时，对应阶模态响应在全频域上减小到可以忽略不计，而其他阶模态响应的在此频率处的峰值也完全消失。

图10 间断附近的模态响应曲线Fig.10 Model response nearby discontinuities

3.2 阻尼对解曲线的影响

分别取ζ1={0.000 1,0.001,0.01}，同时，对于不同的阻尼将频率比的搜索范围扩大为1/β∈(0,120)，得到的解曲线如图11所示。

图11 不同阻尼的解曲线Fig.11 Graphs of solution with different damping

从图11可知，随着阻尼的增大，高阶模态的间断会扩大，直至响应峰值完全消失。需要说明的是：模态响应曲线的峰值并不是和解曲线同时消失的。比如，在图11(b)中，当a=1时，在1/β∈(60,100)区域内看不到不明显的峰值，而此时的前三阶模态响应曲线却如图12所示。

在图12中，1/β∈(60,100)区域内仍能够看到确实存在的峰值这是由于高阶模态阻尼会比低阶模态阻尼大很多，而通过式(45)求解的解曲线的误差也会增大。因而可以说，当解曲线消失时，如若再增大些许阻尼，对应的模态响应峰值也会很快消失。但无论如何，解曲线的消失对于预测频响峰值的消失是有参考价值的。

图12 ζ1=0.001,a=1,1/β∈(0,120)前三阶模态响应曲线Fig.12 First three orders model response (ζ1=0.001,a=1,1/β∈(0,120))

4 结 论

(1)利用第一类拉格朗日方程建立了欧拉-伯努利梁上受单点横向非定常位移约束时的动力学方程，方程为具有时变约束的微分代数方程(DAE)。当非定常约束为谐波函数时，可以得到方程的解析解。使用数值方法验证了解析解的正确性，当约束作用于梁端点时，频响曲线与经实验验证的位移影响函数法相吻合。

(2)利用所得到的解析解，分析了约束工况和力工况的梁的模态响应，研究表明：两者之间的响应具有很大的差别，在约束工况下，在各阶固有频率附近，梁对应阶模态响应不会达到峰值，而其余模达到极小值；当有一个模态响应达到峰值时，其余模态响应也会达到峰值。

(3)利用推导出的求解响应峰值频率比的隐式表达式，可以方便地得到表示非定常约束作用位置与梁动态响应峰值和极小值之间关系的解曲线。对解曲线的分析表明：当非定常约束作用位置从各阶模态函数的零点附近移动到模态函数的零点时，模态响应的峰值和极小值会发生合并，并导致对应阶模态响应在全频域上消失；当阻尼逐渐增大时，高频部分的响应峰值会逐渐消失，各模态的频响峰值和极小值也会消失。

本文所提出的解法及分析研究结论，可用于求解梁的在单点横向非定常约束作用下实际共振频率和稳态响应情况，也为梁在该工况下的设计、分析、测试等问题提供了理论依据。

(A1)

于是

(A2)

其中,

(A3)

由式(A3)可见，该项是仅与β和ζr有关的函数。同样有

(A4)

其中，

(A5)

也是仅与β和ζr有关的函数。

由归一化的模态函数表达式(4)设

(A6)

将(A5)与(A3)代入式(27)中得到

κA=-fl2ω2m2ω32(β,ζr)

(A7)

其中，

ω32(β,ζr)=

(A8)

同样是仅与β和ζr有关的函数。同理，可得

κB=fl2ω2m2ω31(β,ζ1)

(A9)

其中，

ω31(β,ζr)=

(A10)

也是仅与β和ζr有关的函数。

将(A8)与(A10)代入式(20)得到

Ar=ψr(al)fl1.5m0.5ω212,r(β,ζr)

(A11)

其中，

ω212,r(β,ζr)=ω31(β,ζr)ω23,r(β,ζr)-
ω32(β,ζr)ω13,r(β,ζr)

(A12)

是仅与β和ζr有关的函数；同理，有

Br=ψr(al)fl1.5m0.5ω313,r(β,ζr)

(A13)

其中，

ω313,r(β,ζr)=ω31(β,ζr)ω13,r(β,ζr)+
ω32(β,ζr)ω23,r(β,ζr)

(A14)

是仅与β和ζr有关的函数。因此，梁的稳态响应为

(A15)