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

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

(1.天津工業(yè)大學(xué) 機械工程學(xué)院，天津 300387； 2.天津理工大學(xué) 工程訓(xùn)練中心，天津 300384；3.天津理工大學(xué) 機械工程學(xué)院，天津 300384； 4.天津大學(xué) 機械工程學(xué)院，天津 300072)

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

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

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

1 物理建模

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

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

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

(1)

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

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

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

(3)

其中微元段軸向位移為

(4)

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

(5)

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

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

(6)

懸臂梁振動時的彎矩為

(7)

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

梁的彎曲勢能為

(8)

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

(9)

系統(tǒng)的耗散函數(shù)D可表示為

(10)

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

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

(11)

利用Lagrange方程

(12)

(13)

(14)

2 振型函數(shù)和線性基頻

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

(15)

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

(16)

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

(17)

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

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

(18)

假定η(t)=c1cos(β2t)+c2sin(β2t), 其中β2為梁的固有圓頻率,c1，c2為未知的參數(shù)，則式(18)可變?yōu)?/p>

(19)

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

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

(20)

(21)

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

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

(22)

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

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

(23)

式中：Ci為待定系數(shù)；i為模態(tài)的階次，則

(24)

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

(25)

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

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

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

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

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

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

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

γiχicos(Ωt)

(26)

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

(27)

表4 楔形梁前兩階模態(tài)的方程系數(shù)Tab.4 The coefficients of first two modes for wedge beam

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

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

(28)

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

(29)

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

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

(30)

(31)

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

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

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

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

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

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

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

(32)

表8 錐形梁前兩階模態(tài)的方程系數(shù)Tab.8 The coefficients of first two modes for wedge beam

3 主共振響應(yīng)

3.1 多尺度法

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

在主共振情況下，假定

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

(33)

式中，σ為調(diào)諧參數(shù)。將式(33)代入式(26)中可得

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

(34)

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

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

(35)

引入如下微分算子

(36)

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

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

(37)

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

(38)

(39)

(40)

3.2 幅頻響應(yīng)

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

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

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

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

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

4 結(jié) 論

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

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

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