基于Bessel和Meijer-G函數的楔形和錐形懸臂梁振動分析

周坤濤， 楊 濤， 葛 根， 郝淑英， 張琪昌

(1.天津工業大學 機械工程學院，天津 300387； 2.天津理工大學 工程訓練中心，天津 300384；3.天津理工大學 機械工程學院，天津 300384； 4.天津大學 機械工程學院，天津 300072)

變截面梁具有優化質量和強度分布等特點，近年來被廣泛的應用在航空、機械、建筑、橋梁及能量采集器等領域，隨著科學技術的快速發展，變截面梁的設計與加工更趨成熟，為變截面梁的廣泛應用提供了可能。與等截面梁相比，變截面梁的控制方程為四階變系數偏微分方程，通常情況下很難得到解析解[1-2]。

國內外學者對變截面梁進行了大量研究，Silva等[3]將變截面梁假定為厚度不變，寬度線性變化的截面梁，采用四個Meijer-G函數的線性組合表示振型函數，得到了梁的基頻和振型，發現線性問題近似解中可忽略不計的誤差也會對非線性問題的分析結果產生重要影響，進一步強化了線性基頻和振型精確解的重要性。國內學者崔燦等[4]基于歐拉-伯努利梁理論，采用分段傳遞矩陣法對變截面梁的振動進行了研究，該方法的精度依賴于分段數的選??；周坤濤等[5-6]利用超幾何函數和Meijer-G函數的線性組合表示振型函數，對變截面進行自由振動與非線性振動研究，通過數值模擬與試驗方法驗證了其理論。然而，以上研究選取的變截面梁面積和慣性矩均為線性變化。針對寬度不變厚度線性變化的楔形梁以及寬度和厚度均線性變化的錐形梁國內外學者也進行了大量研究，取得了豐碩成果。眾所周知，當寬度不變厚度線性變化時，其面積呈線性變化而慣性矩會呈3次方系數變化；當寬度與厚度均發生線性變化時，其面積會呈2次方系數變化而慣性矩則會呈4次方系數變化，這給解析解的求解帶來很大的困難。針對這一求解難點，Mabie等[7-10]均基于貝塞爾函數理論對楔形或錐形梁進行了研究，雖然其設解的形式不同，最終都得到了具有經典邊界條件的楔形和錐形變截面梁的精確頻率方程及基頻，但是該方法嚴重依賴于方程的形式，方程必須為規范的Bessel方程；Rao[11]采用伽遼金方法得到了變截面梁的基頻；Conway等[12]利用近似多項式表示貝塞爾函數計算錐形梁和楔形梁在簡支、固支和自由邊界條件下的近似基頻；Naguleswaran[13]采用Frobenius級數方法給出了歐拉-伯努利楔形梁和錐形梁的直接解，并以列表的形式給出了精確的基頻；Gaines等[14]考慮剪切變形和轉動效應，利用瑞利-里茲法對錐形和楔形變截面梁前三階固有頻率進行了計算；Lee等[15]利用格林函數進行拉普拉斯變換分析了變截面梁的振動。Ho等[16]采用微分變換法研究了彈性約束邊界條件下變截面梁的振動。Hsu等[17]采用AMDM(Adomain modified decomposition method)法對楔形和錐形變截面梁進行振動分析，得到的線性基頻和振型具有很高的精度，但該方法在四階積分時由于非線性幾何函數會產生奇異性，一般需要進行級數展開，同時求解時需要進行大量的迭代計算，比較耗費機時，其結果的收斂性在很大程度上取決于級數截斷項數的選取。近年來，Lee等[18]對楔形和錐形歐拉-伯努利梁的偏微分方程采用Frobenius冪級數進行設解，提出一種求解楔形和錐形變截面梁的自由振動特性的傳遞矩陣法。Keshmiri等[19]采用改進的數學方法，利用ADM(Adomain decomposition method)法推導了非線性錐形截面歐拉-伯努利梁的特征方程和模態函數，給出了具有不同錐比的指數函數錐形梁和三角函數錐形梁的固有頻率和振型，但該方法在求解過程中依然存在高階泰勒展開，導致計算迭代過程會不斷疊乘，計算仍然耗費機時。李偉等[20]采用微分變換法(differential transform method)推導了圓形變截面梁橫向振動偏微分方程的級數解，并利用有限元軟件ABAQUS驗證了所提方法的精確性。

針對楔形和錐形變截面梁已有的研究成果，本文提出一種無需迭代及近似截斷的Bessel函數與Meijer-G函數線性組合的振型函數，該方法不依賴于楔形和錐形變截面梁的彎曲振動的運動方程是否為標準的Bessel形式，可直接利用Mathematica等計算軟件進行求解，通過給定相應的邊界條件即可快速求解歐拉-伯努利楔形和錐形梁的線性基頻和振型；隨后將本文的振型代入非線性振動控制方程中，考慮梁幾何非線性和慣性非線性的影響，利用多尺度法，研究楔形和錐形梁在主共振下的非線性幅頻響應，Hsu等研究中的計算結果驗證了本文提出的求解振型函數方法的精確性和可靠性。

1 物理建模

如圖1所示的歐拉-伯努利梁為厚度和寬度沿長度方向均逐漸變窄的懸臂梁，彈性模量為E，密度為ρ，長度為L，固定端寬度為b0，厚度為h0，自由端寬度為bl，厚度為hl。建立圖1(a)所示直角坐標系，x軸位于梁的中性軸，y軸沿梁厚度方向，z軸沿梁寬度方向，s軸為沿梁長度方向固定在中性軸上的弧坐標。

圖1 歐拉-伯努利梁理論Fig.1 The Euler-Bernoulli beam theory

定義梁寬度和厚度方向上的截面變化系數分別為：βb=1-bl/b0，βh=1-hl/h0，則截面寬度b(s)和厚度h(s)可以表示為

(1)

由式(1)可以得到坐標s處梁的橫截面積A(s)和截面慣性矩I(s)

假定梁固定端受橫向簡諧位移激勵h(τ)作用產生非線性彎曲振動，忽略梁的重力效應和轉動效應。圖1(b)為坐標s處微段ds的變形圖，該微元段變形包含沿x軸水平位移u(s,τ)和沿y軸的橫向位移w(s,τ)，圖中微元段的位移可表示為

d(s,τ)=u(s,τ)i+[w(s,τ)+h(τ)]j

(3)

其中微元段軸向位移為

(4)

式(4),ξ為虛擬符號。對式(3)求一階導數，并考慮式(4)，即可得到微元段的速度vP。

(5)

式中，(·)為對時間τ求偏導。

當懸臂梁在外激作用下振動時，其動能T為

(6)

懸臂梁振動時的彎矩為

(7)

式中,(′)為對s求偏導。

梁的彎曲勢能為

(8)

將式(7)代入式(8)中，略去高階小量，可得

(9)

系統的耗散函數D可表示為

(10)

假定第i階梁的位移表示為

wi(s,τ)=φi(s)qi(τ)

(11)

利用Lagrange方程

(12)

(13)

(14)

2 振型函數和線性基頻

為確定微分方程式(14)中的系數項，需要求得滿足正交條件和邊界條件的模態函數φi(?)。為了得到梁的頻率和模態函數，首先要求解梁的線性特征值問題。由Euler-Bernoulli理論，該非均勻彈性梁橫向自由振動的偏微分方程

(15)

考慮懸臂梁的邊界條件，固定端約束處時撓度與轉角分別為零，即s=0可得

(16)

在自由端彎矩與剪力分別為零，即s=L可得

(17)

采用相同的無量綱，將式(11)代入式(15)中，可得

[(1-βb?)(1-βh?)3φ″(?)]″η=0

(18)

假定η(t)=c1cos(β2t)+c2sin(β2t), 其中β2為梁的固有圓頻率,c1，c2為未知的參數，則式(18)可變為

(19)

無量綱邊界條件可以重新表示為

φ(0)=0,φ′(0)=0

(20)

(21)

2.1 Euler-Bernoulli楔形梁(βb=0,βh=α)

考慮梁寬度不變，厚度沿x方向線性漸細變化的楔形梁，則式(19)可簡化為

(22)

為求解式(22)的解，將振型函數φ(?)直接表示成第一類Bessel函數與Meijer-G函數線性組合的形式

φi(?)=C1iφ1i+C2iφ2i+C3iφ3i+C4iφ4i

(23)

式中：Ci為待定系數；i為模態的階次，則

(24)

為確定式(23)中的待定系數Ci和βi，需考慮式(20)和式(21)可得

(25)

表1 楔形梁前三階模態的系數βi和Tab.1 The βiand coefficients of the first three modes for wedge beam

圖2 楔形梁前三階振型Fig.2 First three modal shape for wedge beam

為驗證本文方法的正確性，將計算所得的結果與文獻結果進行比對，通過比對表2數值可以看出，本文計算結果與已有文獻結果高度吻合，充分說明本文方法具有很好的計算精度。

表2 楔形懸臂梁前三階無量綱固有圓頻率Tab.2 Non-dimensional first three natural frequencies for a linearly tapered cone cantilever beam

為進一步說明本文所提方法在計算時間和收斂性方面的優勢，選取Windows10操作系統，其中CPU為i5-7300HQ,內存為8 G的筆記本電腦，采用Mathematica11.2軟件編制計算程序，以截面系數α=0.5的楔形梁為例，比對了本文方法和Hsu等方法的計算時間，如表3所示。

表3 截面系數α=0.5的楔形梁前兩階無量綱固有頻率計算時間Tab.3 Calculation time for first two natural frequencies of tapered wedge beams with ratio α=0.5

從表7中可以看出，本文方法能直接得到精確的計算結果，無需考慮展開項數，而AMDM法的計算精度則嚴重依賴于級數項數n的選取，由Hsu等的研究可知，一般n至少超過27時才能得到較精確的結果。從表7中可以發現，求解一階固有頻率時，本文方法只需25.86 s，Hsu等的研究則需要1 375.25 s，二階固有頻率本文方法計算時間為32.02 s，而文獻[17]則需要1 339.28 s。通過比較兩種方法的計算時間，證明本文方法具有計算時間短，收斂快等特點。本文方法在求解截面系數α≤0.5的梁第三階固有頻率時會出現計算誤差或者求解不成功的現象，因此沒有給出第三階固有頻率的具體結果。但對于α>0.5的楔形梁前三階固有頻率的求解，該方法具有求解快，精度高，適用性好等優勢。

γiχicos(Ωt)

(26)

式中:μi為無量綱的阻尼系數；ωi為無量綱的固有頻率;α1i為無量綱彎曲非線性項系數；α2i為無量綱慣性非線性項系數；γi為無量綱外部激勵系數，具體系數為

(27)

表4 楔形梁前兩階模態的方程系數Tab.4 The coefficients of first two modes for wedge beam

2.2 Euler-Bernoulli錐形梁(βb=βh=α)

考慮梁寬度和厚度沿x方向按相同的截面系數線性漸細變化的楔形梁，則式(19)可簡化為

(28)

采用與上文相同的方法，將振型函數φ(?)表示成第二類Bessel函數與Meijer-G函數線性組合的形式，即

(29)

考慮式(20)和式(21)，其邊界條件可表示為

φ(0)=0,φ′(0)=0

(30)

(31)

將式(29)代入式(23)后，考慮式(30)和式(31)，按照2.1節中的方法即可得到振型函數系數Ci和βi，如表5所示。

表5 錐形梁前三階模態的系數βi和Tab.5 The βi and coefficients of the first three modes for cone beam

為驗證本文理論的計算精度，選取不同截面系數的錐形懸臂梁計算了前三階無量綱固有圓頻率，如表6所示，從表中可以看出，一階固有頻率隨著截面系數的增加逐漸變大，而二階和三階逐漸變小。

表6 錐形變截面懸臂梁前三階無量綱固有圓頻率Tab.6 Non-dimensional first three natural frequencies for a linearly tapered cone cantilever beam

同理，為說明本文方法在求解時計算時間和收斂性方面的優勢，選取α=0.7的錐形變截面梁為例，比對了該方法和Hsu等的方法前二階固有頻率的計算時間，如表7所示。

表7 截面系數α=0.7的錐形梁前兩階固有頻率計算時間Tab.7 Calculation time for first two natural frequencies of tapered cone beams with ratio α=0.7

從表7中可以看出，與楔形截面梁計算規律一致，本文方法在處理錐形截面梁固有頻率時同樣具有計算時間短，收斂快等特點，而AMDM法計算精度收斂性則依然嚴重依賴級數項數的選取，但時間會成倍增加。

(32)

表8 錐形梁前兩階模態的方程系數Tab.8 The coefficients of first two modes for wedge beam

3 主共振響應

3.1 多尺度法

為探討變截面梁在主共振下的動態響應，假定外激頻率接近模態頻率(Ω≈ω)且不考慮內共振的情況，采用多尺度法計算變截面懸臂梁的幅頻響應。

在主共振情況下，假定

cos(Ωit)=cos[(ωi+σ)t]

(33)

式中，σ為調諧參數。將式(33)代入式(26)中可得

γiχicos[(ωi+σ)t]

(34)

將時間變量展開成多個不同尺度的時間變量

Tn=εnt,n=0,1,2,…

(35)

引入如下微分算子

(36)

將模態坐標ηn(t)展開成如下形式

ηn(t;ε)=η0(T0,T1)+εη1(T0,T1)+ο(ε2)

(37)

將式(37)代入式(34)，并考慮非線性項，阻尼項，力項均為同階小量，令ε0,ε1的系數為零可得

(38)

(39)

(40)

3.2 幅頻響應

以βb=0，βh=α=0.5楔形梁為例，假設無量綱外激加速度χ=0.019 4，阻尼比μ=0.05，將表3中參數代入式(40)，可得前二階主共振幅頻響應曲線，從圖3中可以看出，一階幅頻響應曲線向右偏，呈現硬特性；二階幅頻響應曲線向左偏，呈現軟特性。

圖3 βh=α=0.5楔形梁在外激作用下的幅頻響應曲線Fig.3 Amplitude-frequency response curve of wedge beam under external excitation for βh=α=0.5

為驗證本文幅頻響應曲線的正確性，選取Hsu等所用的AMDM法得到的振型函數，將其代入式(27)，得到相應的系數，隨后代入式(40)同樣可以得到幅頻響應曲線，從圖中可以看出，兩種方法所得到的幅頻響應曲線高度吻合。眾所周知，線性基頻和振型函數微小的誤差，將會對非線性分析產生很大的影響，這也充分說明本文方法得到的振型函數與AMDM法得到的振型函數高度吻合。

同理，以βb=βh=α=0.7的錐形梁為例，采用同樣的外激加速度作用，得到了前兩階主共振幅頻響應曲線，從圖4中可以看出，一階幅頻響應曲線兩種方法吻合很好，二階幅頻響應曲線兩種方法僅有微小的誤差，其原因可能是AMDM法振型函數的精確性與級數展開及截斷項數的選取有關，而本文理論方法則不需要考慮近似截斷項數。

圖4 βb=βh=α=0.7錐形梁在外激作用下的幅頻響應曲線Fig.4 Amplitude-frequency response curve of cone beam under external excitation for βb=βh=α=0.7

4 結 論

(1)本文提出了一種快速求解楔形和錐形截面懸臂梁線性基頻和模態函數的新方法，該方法不依賴于楔形和錐形截面梁的彎曲振動的運動方程是否為標準的Bessel形式，直接利用Bessel函數和Meijer-G函數線性組合構建新的振型函數，求解過程無需大量迭代運算和級數近似截斷。通過理論計算得到的線性基頻與已有文獻結果高度吻合，驗證了本文方法的正確性。

(2)線性基頻和振型函數微小的誤差，將會對非線性幅頻響應分析產生很大的影響。從本文得到的幅頻響應曲線圖可以看出：采用Bessel函數和Meijer-G函數線性組合構建新的振型函數與AMDM法得到的振型函數一樣，具有很高的精度，驗證了本文振型函數的正確性。

(3)本文方法在處理變截面參數較大的楔形和錐形梁時有更好的適用性，該方法可為楔形和錐形截面梁固有頻率及振型函數的求解提供了新的思路。