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

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

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

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

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

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

1 物理建模

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

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

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

(1)

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

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

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

(3)

其中微元段軸向位移為

(4)

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

(5)

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

當(dāng)懸臂梁在外激作用下振動(dòng)時(shí)，其動(dòng)能T為

(6)

懸臂梁振動(dòng)時(shí)的彎矩為

(7)

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

梁的彎曲勢(shì)能為

(8)

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

(9)

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

(10)

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

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

(11)

利用Lagrange方程

(12)

(13)

(14)

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

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

(15)

考慮懸臂梁的邊界條件，固定端約束處時(shí)撓度與轉(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方向線性漸細(xì)變化的楔形梁，則式(19)可簡(jiǎn)化為

(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

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

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

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

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

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

γiχicos(Ωt)

(26)

式中:μi為無量綱的阻尼系數(shù)；ωi為無量綱的固有頻率;α1i為無量綱彎曲非線性項(xiàng)系數(shù)；α2i為無量綱慣性非線性項(xiàng)系數(shù)；γi為無量綱外部激勵(lì)系數(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ù)線性漸細(xì)變化的楔形梁，則式(19)可簡(jiǎn)化為

(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

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

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

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

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

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

(32)

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

3 主共振響應(yīng)

3.1 多尺度法

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

在主共振情況下，假定

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

(33)

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

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

(34)

將時(shí)間變量展開成多個(gè)不同尺度的時(shí)間變量

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

(35)

引入如下微分算子

(36)

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

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

(37)

將式(37)代入式(34)，并考慮非線性項(xiàng)，阻尼項(xiàng)，力項(xiàng)均為同階小量，令ε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àn)證本文幅頻響應(yīng)曲線的正確性，選取Hsu等所用的AMDM法得到的振型函數(shù)，將其代入式(27)，得到相應(yīng)的系數(shù)，隨后代入式(40)同樣可以得到幅頻響應(yīng)曲線，從圖中可以看出，兩種方法所得到的幅頻響應(yīng)曲線高度吻合。眾所周知，線性基頻和振型函數(shù)微小的誤差，將會(huì)對(duì)非線性分析產(chǎn)生很大的影響，這也充分說明本文方法得到的振型函數(shù)與AMDM法得到的振型函數(shù)高度吻合。

同理，以βb=βh=α=0.7的錐形梁為例，采用同樣的外激加速度作用，得到了前兩階主共振幅頻響應(yīng)曲線，從圖4中可以看出，一階幅頻響應(yīng)曲線兩種方法吻合很好，二階幅頻響應(yīng)曲線兩種方法僅有微小的誤差，其原因可能是AMDM法振型函數(shù)的精確性與級(jí)數(shù)展開及截?cái)囗?xiàng)數(shù)的選取有關(guān)，而本文理論方法則不需要考慮近似截?cái)囗?xiàng)數(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ù)的新方法，該方法不依賴于楔形和錐形截面梁的彎曲振動(dòng)的運(yùn)動(dòng)方程是否為標(biāo)準(zhǔn)的Bessel形式，直接利用Bessel函數(shù)和Meijer-G函數(shù)線性組合構(gòu)建新的振型函數(shù)，求解過程無需大量迭代運(yùn)算和級(jí)數(shù)近似截?cái)?。通過理論計(jì)算得到的線性基頻與已有文獻(xiàn)結(jié)果高度吻合，驗(yàn)證了本文方法的正確性。

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

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