陳美霞, 姚仕輝, 謝 坤
(華中科技大學(xué) 船舶與海洋工程學(xué)院,武漢 430074)
圓板作為常見的結(jié)構(gòu)單元,在工程上有著廣泛的應(yīng)用。在很多情況下,圓板是與流體接觸的,如微量泵、蝶閥的圓盤、整體式核反應(yīng)堆的冷卻劑等。流體的存在會(huì)顯著降低圓板的固有頻率,從而改變結(jié)構(gòu)的動(dòng)力響應(yīng)特性,因此對(duì)圓板和流體的流-固耦合問題進(jìn)行研究具有重要意義。
Kwak[1]研究了漂浮在無限流域上圓板的振動(dòng)特性,通過引入無量綱附加質(zhì)量因子(NAVMI)來評(píng)價(jià)流體對(duì)于圓板的作用,近似求解了圓板的固有頻率。后來,Kwak[2]應(yīng)用Hankel變換計(jì)算了嵌在剛性障板中圓板的水動(dòng)力特性,通過與文獻(xiàn)[1]進(jìn)行比較,發(fā)現(xiàn)嵌在剛性障板中的圓板所受到的流體附加質(zhì)量作用比漂浮在自由液面上的圓板更大。Kwak等[3]求解了漂浮在深度有限的流體上圓板的振動(dòng)特性,并分析了深度改變對(duì)于圓板固有頻率的影響。Kwak等通過NAVMI法求解圓板的固有頻率時(shí),忽略了流體對(duì)圓板振型的影響,因此計(jì)算結(jié)果有一定誤差。Amabili等[4]在考慮流體對(duì)結(jié)構(gòu)振型影響的基礎(chǔ)上,通過Rayleigh-Ritz法建立控制方程,更精確地求解了結(jié)構(gòu)的固有頻率,并與NAVMI法進(jìn)行對(duì)比,結(jié)果表明,采用NAVMI法求解低階頻率時(shí)誤差較小。Bauer[5]研究了彈性薄膜封蓋或彈性薄板封蓋與容器內(nèi)流體的耦合振動(dòng),結(jié)果表明,封蓋與流體耦合振動(dòng)的固有頻率會(huì)大于沒有封蓋時(shí)流體的晃動(dòng)頻率。Bauer等[6]在Bauer的研究基礎(chǔ)上,考慮了流體粘性的影響。Amabili[7]研究了漂浮在容器內(nèi)流體表面的圓板的振動(dòng)特性,分析了自由液面表面波對(duì)于圓板振動(dòng)的影響。在文獻(xiàn)[7]的基礎(chǔ)上,Yousefzadeh等[8]分析了漂浮在自由液面上的功能梯度圓板的振動(dòng)特性。Cheung等[9]研究了作為開口容器底部結(jié)構(gòu)的彈性圓板的水動(dòng)力特性,利用疊加原理將自由液面邊界條件分離出來,分析了流體自由表面波對(duì)圓板振動(dòng)的影響。Jeong等[10]研究了充滿可壓縮流體的封閉容器中圓形隔板的振動(dòng)特性,分析了流體壓縮性對(duì)于圓板振動(dòng)特性的影響。Askari等[11]對(duì)浸沒在開口容器中的圓板的振動(dòng)特性進(jìn)行分析,并考慮了液體晃動(dòng)的影響。Eftekhari[12]引入微分求積法(DQM)求解了圓板的流-固耦合振動(dòng)問題。Tariverdilo等[13]分別用Fourier-Bessel級(jí)數(shù)法和變分法分析了圓板的流-固耦合問題,并對(duì)比了不同方法的計(jì)算結(jié)果,但未涉及軸對(duì)稱模態(tài)的計(jì)算。Jeong[14]計(jì)算了完全相同的兩塊圓板與流體的耦合振動(dòng),分別求解了耦合系統(tǒng)的同相和反相模態(tài)。在求解軸對(duì)稱模態(tài)時(shí),Jeong對(duì)附加質(zhì)量矩陣的元素進(jìn)行了修正,但反相軸對(duì)稱模態(tài)結(jié)果的誤差仍相對(duì)較大。
基于文獻(xiàn)[13-14],本文結(jié)合Fourier-Bessel級(jí)數(shù)展開法和Galerkin法建立控制方程,對(duì)與流體接觸的邊界彈性約束圓板的固有振動(dòng)特性進(jìn)行了研究。利用0階貝塞爾函數(shù)的性質(zhì),添加附加約束方程,對(duì)軸對(duì)稱模態(tài)進(jìn)行求解。在此基礎(chǔ)上分析了彈簧剛度和流體深度對(duì)于圓板自由振動(dòng)的影響。
如圖1所示,圓柱形容器內(nèi)充滿不可壓縮的無黏性流體,容器底面和柱面為剛性壁面,頂面覆有彈性圓板。流體密度為ρ0,深度為d;圓板半徑為a,板厚為h,圓板邊緣設(shè)有均勻分布的位移彈簧和轉(zhuǎn)角彈簧,彈簧的分布剛度分別為K和Kφ。
假設(shè)圓板滿足基爾霍夫薄板理論,則系統(tǒng)的動(dòng)力平衡方程可表示為
(1)
式中:w,D,ρ分別為板的撓度、彎曲剛度和密度,P為流體壓力,▽4為柱坐標(biāo)系下的雙調(diào)和算子,其表達(dá)式為
(2)
根據(jù)模態(tài)正交性理論,圓板位移可以展開為真空中振型的級(jí)數(shù):
(3)
圖1 與流體接觸的邊界彈性約束圓板Fig.1 Elastically restrained circular plate in contact with fluid
式中:Wnm為圓板在真空中的振型,qnm為廣義坐標(biāo),n和m分別為節(jié)徑數(shù)和節(jié)圓數(shù)。Wnm滿足式(1)對(duì)應(yīng)的齊次微分方程:
(4)
根據(jù)文獻(xiàn)[15]可知,圓板在真空中的振型為:
Wnm(r,θ)=[Jn(λnmr)+
CnmIn(λnmr)]cos(nθ)
(5)
式中:Jn和In分別為n階第一類貝塞爾函數(shù)和第一類修正貝塞爾函數(shù)。Cnm為與振型有關(guān)的系數(shù),λnm為頻率系數(shù),滿足
(6)
式中:ωnm表示真空中圓板的固有頻率。
當(dāng)圓板材料和尺寸一定時(shí),振型系數(shù)Cnm和頻率系數(shù)λnm取決于彈簧剛度K和Kφ,系數(shù)的求解步驟見附錄A。由此可知,當(dāng)圓板邊界條件確定時(shí),真空中的振型Wnm為已知函數(shù),只要求解出廣義坐標(biāo)qnm,即可根據(jù)式(3)得到圓板在流體中的振型。
假設(shè)流體作無旋運(yùn)動(dòng),則可用速度勢(shì)函數(shù)Φ(r,θ,x,t)描述流體運(yùn)動(dòng)。結(jié)合1.1節(jié)中流體不可壓縮且無粘性的假設(shè),Φ(r,θ,x,t)滿足拉普拉斯方程:
▽2Φ=0
(7)
將速度勢(shì)函數(shù)的時(shí)間項(xiàng)和空間項(xiàng)分離,可得[14]
Φ(r,θ,x,t)=iωφ(r,θ,x)eiωt
(8)
式中:φ為空間速度勢(shì)函數(shù),且同樣滿足拉普拉斯方程。應(yīng)用分離變量法可得,對(duì)于某一個(gè)確定的節(jié)徑數(shù)n,φ的一般解為(參考附錄B):
φ(r,θ,x)=(anx+bn)δn0cosnθ+
Fnscosh(βnsx)]cosnθ
(9)
式中:δn0是Kronecker delta函數(shù),an、bn、Ens、Fns、βns為待求系數(shù),且βns>0。
由于容器柱面和底面為剛性壁面,φ需要滿足如下邊界條件:
(10)
(11)
將式(9)代入(10),可以得到關(guān)于βns的方程:
(12)
同理,將式(9)代入(11),可得
a0=0
(13)
Ens=0
(14)
因此,式(9)可以簡(jiǎn)化為
φ(r,θ,x)=cosnθ[b0δn0+
(15)
在流-固耦合面,需要滿足圓板和流體的速度連續(xù)性條件:
(16)
將式(3)、(8)和(15)代入式(16)可得
(17)
其中Rnm(r)=Jn(λnmr)+CnmIn(λnmr)為干圓板振型的徑向分量。根據(jù)貝塞爾方程本征函數(shù)的正交性[16],在方程(17)左右兩邊乘以rJn(βnsr),并在[0,a]積分可得
(18)
當(dāng)n=0時(shí),常數(shù)函數(shù)也是貝塞爾方程的本征函數(shù),且對(duì)于?β0s>0,常數(shù)函數(shù)與J0(β0sr)在[0,a]帶權(quán)r正交。在方程(17)兩邊同時(shí)乘以r并在[0,a]積分,可得n=0時(shí)還需要滿足的附加約束方程:
(19)
將式(18)代入(15),速度勢(shì)函數(shù)φ可以進(jìn)一步表示為
φ(r,θ,x)=cosnθ[b0δn0+
(20)
其中
(21)
假設(shè)圓板作小振幅振動(dòng),應(yīng)用線性伯努利方程可得作用于圓板的流體壓力:
(22)
此時(shí)對(duì)于確定的n,平衡方程可以表示為
(23)
應(yīng)用Galerkin法,在方程(23)兩邊乘以Wnk(r,θ)并在圓板范圍積分,結(jié)合真空中振型的正交性[15]可得代數(shù)方程:
(24)
分別取m和s的截?cái)鄶?shù)為M和S,可以將式(24)寫成矩陣形式。當(dāng)n>0時(shí),式(24)可表示為
[ρhPn-ω2(ρhZn+ρ0Gn)]qn=0
(25)
其中qn={qn0,qn1,…,qnM}T為系數(shù)向量,Pn,Zn,Gn為M+1階矩陣,各矩陣元素分別為
(26)
(27)
k,m=0,1,2,…,M
(28)
當(dāng)n=0時(shí),除系數(shù)向量q0外b0也是未知量,為求得b0和q0,代數(shù)方程組(24)需要和附加約束方程(19)聯(lián)立求解,表示成矩陣形式分別為
[ρhP0-ω2(ρhZ0+ρ0G0)]q0-ρ0ω2c0b0=0
(29)
l0q0=0
(30)
式中:P0,Z0,G0和q0與n>0的情況有著相同的表達(dá)式。c0和l0分別為M+1維的列向量和行向量,其元素分別為
(31)
(32)
此時(shí)n=0時(shí)的控制方程可以表示為如下矩陣形式:
(33)
式(25)和(33)給出了矩陣形式的控制方程。要使方程存在非平凡解,必須讓系數(shù)矩陣的行列式為零,由此求解出的ω即為與流體接觸的圓板的固有頻率;求出對(duì)應(yīng)的非零解即可得廣義坐標(biāo)qnm,代入式(3)可得圓板在流體中的振型。
為驗(yàn)證本文解析方法的正確性,將基于本文解析方法的計(jì)算結(jié)果與基于有限元法的數(shù)值結(jié)果進(jìn)行對(duì)比。計(jì)算時(shí)采用如下參數(shù):圓板半徑a=0.2 m,板厚h=0.002 m,密度ρ=2 700 kg/m3,楊氏模量E=69 GPa,泊松比ν=0.3;流體密度ρ0=1 000 kg/m3,深度d=0.02 m;位移彈簧的分布剛度K=103N/m2,轉(zhuǎn)角彈簧的分布剛度Kφ=103N/rad。本文方法采用Matlab實(shí)現(xiàn),截?cái)鄶?shù)M和S均取40。數(shù)值結(jié)果采用有限元軟件Ansys得到,其中圓板采用Shell63單元,流體采用Fluid30單元,網(wǎng)格尺寸取Δ=0.004 m。
對(duì)截?cái)鄶?shù)M和S的收斂性分析如表1所示。由表1可知,在保留兩位小數(shù)時(shí),截?cái)鄶?shù)取40和50的計(jì)算結(jié)果是一致的,而截?cái)鄶?shù)取30時(shí)部分高階頻率尚未收斂,故取截?cái)鄶?shù)M=S=40是合理的。
表1選取不同截?cái)鄶?shù)M和S求得的圓板固有頻率
Tab.1Naturalfrequenciesofthecircularplateobtainedbyusingdifferenttruncationnumber
nm固有頻率/HzM=S=30M=S=40M=S=5001118.63118.63118.632646.26646.26646.2631 697.141 697.141 697.141011.6111.6011.601305.95305.94305.9421 085.351 085.351 085.3532 409.992 409.992 409.992058.1958.1958.191577.72577.72577.7221 618.651 618.651 618.6533 220.713 220.703 220.7030153.58153.58153.581932.13932.13932.1322 242.582 242.572 242.5734 126.834 126.824 126.82
為驗(yàn)證有限元網(wǎng)格的收斂性,分別對(duì)Δ1=0.008 m,Δ2=0.004 m,Δ3=0.002 m三種網(wǎng)格尺寸的有限元模型進(jìn)行分析,計(jì)算結(jié)果如表2所示。表中f1、f2、f3為不同模型的固有頻率,e1、e2為相對(duì)誤差,滿足e1=(f1-f2)/f2,e2=(f2-f3)/f3。由于有限元網(wǎng)格收斂較慢,在綜合考慮計(jì)算速度和精度的條件下,本文以相對(duì)誤差<1%作為有限元網(wǎng)格收斂的評(píng)判標(biāo)準(zhǔn),結(jié)合表2的結(jié)果,取網(wǎng)格尺寸Δ=Δ2=0.004 m是合理的。
表2采用不同網(wǎng)格尺寸的有限元模型得到的圓板固有頻率
Tab.2NaturalfrequenciesofthecircularplateobtainedbyFEMusingdifferentmeshsize
nmΔ1=0.008 mΔ2=0.004 mΔ3=0.002 mf1/Hze1/%f2/Hze2/%f3/Hz01119.230.37118.790.11118.662653.790.88648.100.26646.4131 728.601.441 704.090.421 696.911011.630.1511.610.0411.611308.280.57306.530.17306.0221 101.501.141 089.110.341 085.4532 463.001.722 421.270.512 408.992058.420.2958.250.0958.201583.250.72579.080.22577.8121 644.901.261 624.410.381 618.2833 293.901.803 235.520.553 217.8830154.720.55153.870.15153.631944.291.00934.960.30932.1822 288.701.612 252.420.482 241.7034 244.202.264 150.470.684 122.31注:e1=(f1-f2)/f2,e2=(f2-f3)/f3
分別采用本文方法和有限元法求得的圓板固有頻率結(jié)果如表3所示。從表中可知,對(duì)于n,m=0,1,2,3的模態(tài),本文方法與有限元法計(jì)算結(jié)果相對(duì)誤差在1%
表3本文方法與有限元法求得的固有頻率結(jié)果對(duì)比
Tab.3ComparisonofthenaturalfrequenciesobtainedbypresentmethodandFEM
nm有限元法/Hz本文方法/Hz相對(duì)誤差/%01118.79118.63-0.132648.10646.26-0.2831 704.091 697.14-0.411011.6111.60-0.051306.53305.94-0.1921 089.111 085.35-0.3532 421.272 409.99-0.472058.2558.19-0.111579.08577.72-0.2421 624.411 618.65-0.3533 235.523 220.70-0.4630153.87153.58-0.191934.96932.13-0.3022 252.422 242.57-0.4434 150.474 126.82-0.57
以內(nèi),驗(yàn)證了本文方法的正確性。此外,由表3還可以看出隨著節(jié)徑數(shù)n和節(jié)圓數(shù)m的增加,圓板的固有頻率變大。需要注意的是,本文假定流體不可壓縮,因此流體體積保持不變,故不存在n=m=0的模態(tài)。
2.2.1 彈簧剛度對(duì)圓板自由振動(dòng)特性的影響
(a) 自由(K=0, Kφ=0)
(c) 導(dǎo)向(K=0, Kφ=∞)
(d) 固支(K=∞, Kφ=∞)圖2 不同邊界條件下圓板的無量綱頻率Fig.2 Nondimentional frequencies of circular plates with different boundary conditions
流體不僅會(huì)降低圓板的固有頻率還會(huì)改變圓板的振型,圖3給出了圓板在真空中和流體中的歸一化振型的對(duì)比。由圖3看出,簡(jiǎn)支和固支圓板的振型受到流體的影響相對(duì)較大,自由邊界圓板的振型受到流體影響很小,導(dǎo)向邊界圓板在真空中和流體中的振型幾乎完全一致。
(a) n=0,m=2自由
(b) n=1,m=2自由
(c) n=0,m=2簡(jiǎn)支
(d) n=1,m=2簡(jiǎn)支
(e) n=0,m=2導(dǎo)向
(f) n=1,m=2導(dǎo)向
(g) n=0,m=2固支
(h) n=1,m=2固支圖3 圓板在真空中和流體中的歸一化振型Fig.3 Normalized mode shapes of the circular plate in vacuo and in fluid
在求解固有頻率時(shí),若忽略流體對(duì)振型的影響,采用圓板在真空中的振型代替圓板在流體中的振型進(jìn)行求解,即假定w(r,θ)=Wnm(r,θ),式(23)可簡(jiǎn)化為
(34)
當(dāng)n>0時(shí),式(34)兩邊乘以Wnm在[0,a]積分可以直接求出固有頻率ω,相對(duì)于求解矩陣方程(25)的廣義特征值問題,將大為減少計(jì)算時(shí)間。
而對(duì)于n=0,式(19)簡(jiǎn)化為
(35)
此式無法成立,因此在求解軸對(duì)稱模態(tài)(n=0的模態(tài))時(shí),這種簡(jiǎn)化方法是不適用的。
利用式(34)解出的圓板固有頻率如表4所示,表中括號(hào)中結(jié)果為通過求解式(25)的廣義特征值問題得到的固有頻率。由表4可知,對(duì)于自由和導(dǎo)向邊界圓板,忽略流體對(duì)振型的影響對(duì)圓板固有頻率的計(jì)算結(jié)果影響較小,誤差在5%以內(nèi);而對(duì)于簡(jiǎn)支和固支圓板,則會(huì)產(chǎn)生較大的誤差。由圖3的結(jié)果可知,這是因?yàn)樽杂珊蛯?dǎo)向邊界圓板振型受到流體影響更小,采用近似方法計(jì)算時(shí)假定振型(真空中的振型)與真實(shí)振型(流體中的振型)很接近,因此誤差較小。從表4還可以看出,對(duì)于以上四種邊界條件,當(dāng)周向波數(shù)n增加時(shí),近似解法的誤差會(huì)降低。由此可知,當(dāng)邊界條件為自由、導(dǎo)向或者當(dāng)所求模態(tài)周向波數(shù)n較大時(shí),可以忽略流體對(duì)振型的影響。
表4 忽略流體對(duì)振型影響求得的圓板固有頻率Tab.4 Natural frequencies of the circular plate obtained by neglecting the fluid influences on mode shapes
2.2.2 流體深度對(duì)圓板自由振動(dòng)特性的影響
出現(xiàn)上述結(jié)果是由于底部剛性壁面會(huì)限制流體的垂向運(yùn)動(dòng),且距離剛性底面越近,流體作垂向運(yùn)動(dòng)越困難。而圓板受到流體的附加質(zhì)量取決于流-固耦合層附近流體的運(yùn)動(dòng)情況,深度d越小,流-固耦合層附近的流體距離剛性底面越接近,運(yùn)動(dòng)越困難,流體隨圓板振動(dòng)時(shí)需要圓板提供更大的作用力,反過來說圓板受到流體的作用力也就越大,附加質(zhì)量也就越大。
(a) n=0
(b) n=1
(c) n=2
(d) n=3圖4 流體深度對(duì)無量綱頻率的影響Fig.4 The influence of fluid depth on nondimentional frequencies
(1) 本文采用Fourier-Bessel級(jí)數(shù)法和Galerkin法,研究了與流體接觸的邊界彈性約束圓板的自由振動(dòng)特性。針對(duì)軸對(duì)稱模態(tài),利用0階貝塞爾函數(shù)的性質(zhì),添加附加約束方程進(jìn)行求解。本文方法與數(shù)值法結(jié)果吻合良好,驗(yàn)證了本文計(jì)算方法的正確性。
(2) 通過改變彈簧剛度研究了幾種常見邊界條件下圓板的流-固耦合振動(dòng)特性。結(jié)果表明,流體作用會(huì)降低圓板的固有頻率,且對(duì)低階模態(tài)影響更大,此結(jié)論對(duì)于不同邊界條件的圓板均成立。此外,流體還會(huì)影響圓板的振型,其中簡(jiǎn)支和固支邊界圓板的振型受流體影響較大,而自由及導(dǎo)向邊界圓板的振型受到流體影響較小,針對(duì)后兩種邊界條件,可以用圓板在真空中的振型代替圓板在流體中的振型簡(jiǎn)化計(jì)算,且這種簡(jiǎn)化具有較高的精度。
(3) 在較低的深度范圍內(nèi),增大流體深度會(huì)顯著降低流體對(duì)于圓板的附加質(zhì)量作用;當(dāng)深度大于1.5倍圓板半徑時(shí),流體深度的改變對(duì)于圓板自由振動(dòng)的影響可以忽略。
(36)
(37)
式中:Mr、Nr和Mrθ分別代表單位長(zhǎng)度的彎矩、剪力和扭矩,這三種分布內(nèi)力的表達(dá)式為:
(38)
(39)
(40)
式中:ν為圓板的泊松比。將式(38)~(40)和圓板位移表達(dá)式(3)代入式(36)和(37)可得:
h1(λnm)+Cnmh2(λnm)=0
(41)
y1(λnm)+Cnmy2(λnm)=0
(42)
式中:h1(λnm),h2(λnm),y1(λnm)和y2(λnm)為關(guān)于λnm的函數(shù),其表達(dá)式為:
(43)
(44)
(45)
(46)
聯(lián)立式(41)和(42)并消去Cnm,可得真空中圓板的頻率方程:
h1(λnm)y2(λnm)-h2(λnm)y1(λnm)=0
(47)
由式(47)可解得頻率系數(shù)λnm, 將λnm代入式(41)即可得到振型系數(shù)Cnm:
(48)
附錄B流場(chǎng)拉普拉斯方程求解
對(duì)空間速度勢(shì)函數(shù)φ(r,θ,x)采用分離變量法:
φ(r,θ,x)=R(r)Θ(θ)X(x)
(49)
將式(49)代入拉普拉斯方程,并在方程兩邊同時(shí)除以R(r)×Θ(θ)X(x),可得
(50)
由于r,θ,x為相互獨(dú)立的變量,式(50)成立的條件為方程兩邊等于常數(shù):
(51)
其中μ為常數(shù)??傻肵(x)滿足微分方程:
X″+μX=0
(52)
同時(shí),由式(51)可得
(53)
其中λ為常數(shù)??紤]到φ(r,θ,x)在周向以2π為周期,Θ(θ)滿足定解條件:
(54)
由方程(54)可求得常數(shù)λ和函數(shù)Θ(θ)分別為
(55)
在求解自由振動(dòng)特性時(shí),不妨以振型對(duì)稱軸作為r軸,則A2=0,可得
Θ(θ)=A1cosnθ
(56)
將式(56)代入式(53)可知R(r)需要滿足微分方程:
r2R″+rR′-(μr2+n2)R=0
(57)
針對(duì)μ的取值,對(duì)以下幾種情況進(jìn)行分類討論:
1)μ<0,方程(52)的通解為
(58)
此時(shí)方程(57)為n階貝塞爾方程,通解為
(59)
φ(r,θ,x)=Jn(βnsr)cosnθ[Enssinh(βnsx)+
Fnscosh(βnsx)]
(60)
式中:βns>0。
2)μ>0,方程(57)為n階虛宗量貝塞爾方程,通解為
(61)
由R(0)為有限值知D2=0,同時(shí)考慮到柱面剛性壁面條件R′(a)=0得D1=0。因此在μ>0時(shí),φ(r,θ,x)不存在非零解。
3)μ=0,方程(52)的通解為
X(x)=E1x+E2
(62)
此時(shí)方程(57)為歐拉型常微分方程,其通解為
(63)
由R(0)為有限值知F2=G2=0,同時(shí)考慮到R′(a)=0可得F1=0。因此式μ=0時(shí),速度勢(shì)函數(shù)通解為
φ(r,θ,x)=cosnθ(anx+bn)δn0
(64)
綜上可知,對(duì)于確定的n,速度勢(shì)函數(shù)的一般解為
φ(r,θ,x)=cosnθ{(anx+bn)δn0+
(65)