受單點橫向非定常約束梁的響應分析

丁維高，謝 進

(西南交通大學 機械工程學院，成都 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)