具有循環(huán)對稱特點(diǎn)的耦合結(jié)構(gòu)振動特性分析

藺彥虎,趙寧,彭艷軍

(西北工業(yè)大學(xué)機(jī)電學(xué)院,710072,西安)

?

具有循環(huán)對稱特點(diǎn)的耦合結(jié)構(gòu)振動特性分析

藺彥虎,趙寧,彭艷軍

(西北工業(yè)大學(xué)機(jī)電學(xué)院,710072,西安)

利用一些循環(huán)對稱結(jié)構(gòu)件經(jīng)耦合后仍呈現(xiàn)循環(huán)對稱的特點(diǎn)進(jìn)行耦合結(jié)構(gòu)的振動特性分析,以節(jié)約計算時間和計算資源。應(yīng)用Weierstrass-Mandelbrot分型模型結(jié)合赫茲接觸理論求解了耦合結(jié)構(gòu)接觸面的切、法向剛度,再用商業(yè)有限元軟件導(dǎo)出了各結(jié)構(gòu)件的剛度及質(zhì)量矩陣,并通過擴(kuò)維來組裝耦合扇區(qū)的剛度矩陣;應(yīng)用彈簧元模擬了接觸面的切、法向剛度;對耦合結(jié)構(gòu)施加復(fù)約束條件,用MATLAB開發(fā)出利用結(jié)構(gòu)單扇區(qū)求解耦合結(jié)構(gòu)振動特性的計算流程。將該流程應(yīng)用于航空弧齒錐齒輪-阻尼碗耦合結(jié)構(gòu)振動特性的求解,并與耦合結(jié)構(gòu)的整體有限元模型計算結(jié)果進(jìn)行了對比,結(jié)果表明:開發(fā)的計算流程在保證準(zhǔn)確性的前提下,可以高效率地完成結(jié)構(gòu)振動特性分析,計算誤差控制在1%之內(nèi),運(yùn)算時間節(jié)省了67%以上,同時可以有效節(jié)約存儲資源。另外,通過對航空錐齒輪-附加阻尼碗的振動特性分析得知,該耦合結(jié)構(gòu)增加了以阻尼碗振動為主振型的特征頻率,在后續(xù)的結(jié)構(gòu)減振或頻響計算時應(yīng)引起足夠重視。

耦合結(jié)構(gòu);循環(huán)對稱結(jié)構(gòu);振動;特征值;復(fù)約束

旋轉(zhuǎn)結(jié)構(gòu)通常具有循環(huán)對稱的特點(diǎn)。為防止旋轉(zhuǎn)結(jié)構(gòu)件在工作中出現(xiàn)大振動而導(dǎo)致疲勞破壞,一個有效的做法是在該旋轉(zhuǎn)結(jié)構(gòu)件上安裝附加結(jié)構(gòu)件[1]。附加結(jié)構(gòu)件一般對原結(jié)構(gòu)產(chǎn)生2種作用:①調(diào)頻作用;②應(yīng)用結(jié)合面間的干摩擦提供阻尼效應(yīng)。附加結(jié)構(gòu)件的結(jié)構(gòu)一般較為簡單,也呈現(xiàn)出循環(huán)對稱的特點(diǎn)。因此,一般的結(jié)構(gòu)體-附加結(jié)構(gòu)件耦合結(jié)構(gòu)仍具有循環(huán)對稱特征,如航空篦齒封嚴(yán)結(jié)構(gòu)與阻尼襯套(片)[2],航空發(fā)動機(jī)葉片與減振阻尼楔塊,以及薄壁齒輪與阻尼環(huán)、阻尼碗等。附加結(jié)構(gòu)件有可能對原結(jié)構(gòu)的振動形式產(chǎn)生一定影響,因此,對耦合結(jié)構(gòu)進(jìn)行振動特性分析是附加阻尼件減振的基本工作。

單個循環(huán)對稱結(jié)構(gòu)的振動特性求解技術(shù)[3-5]已經(jīng)較為成熟,大型商業(yè)有限元軟件(如ANSYS、Abaqus等)均將該技術(shù)嵌入其求解程序,用戶只需定義單扇區(qū)模型及其主、從界面即可完成整個結(jié)構(gòu)的振動特性求解[6]。然而,對于呈現(xiàn)循環(huán)對稱特點(diǎn)的耦合結(jié)構(gòu)件的振動特性分析,上述軟件卻無法應(yīng)用單扇區(qū)給出相應(yīng)的解。針對此類問題的相關(guān)文獻(xiàn)目前較為鮮見。文獻(xiàn)[2]分析了裝有減振襯套的航空篦齒封嚴(yán)結(jié)構(gòu)的振動特性,但該分析是應(yīng)用整體有限元模型進(jìn)行求解的。

本文應(yīng)用Weierstrass-Mandelbrot(W-M)分形模型結(jié)合赫茲接觸理論確定了耦合界面的切、法向接觸剛度;應(yīng)用Abaqus導(dǎo)出了結(jié)構(gòu)及附加件扇區(qū)的剛度矩陣和阻尼矩陣,并對結(jié)構(gòu)扇區(qū)的剛度及質(zhì)量矩陣進(jìn)行了擴(kuò)維,將附加結(jié)構(gòu)扇區(qū)的剛度及質(zhì)量矩陣裝入擴(kuò)維矩陣,用彈簧元耦合了扇區(qū)的切、法向剛度;應(yīng)用復(fù)約束法求解了組合結(jié)構(gòu)的振動特性。最后以航空錐齒輪-附加阻尼碗為例進(jìn)行了求解,并通過對比求解齒輪-阻尼碗整體模型驗(yàn)證了本文提出的求解流程的有效性及對求解資源的節(jié)約。

1 耦合剛度計算

1.1 W-M分形模型

分形幾何可以描述不同長度范圍的表面幾何粗糙特征,從而形成一種描述粗糙度的方法[7-8]。表面粗糙度可用W-M分形函數(shù)[7-9]來表示

(1)

式中:z為表面形貌高度;x為橫向距離;G為分形粗糙度參數(shù),該參數(shù)獨(dú)立于頻率密度,表示對表面高度的縮放;D為分形維數(shù)(11),控制表面形貌的頻率密度,文獻(xiàn)[8]推薦取γ=1.5;L為幾何形體表面長度。圖1所示為不同分形維數(shù)的物體表面形貌特征。

圖1 不同分形維數(shù)的分形表面輪廓

宏觀平坦的2個平面接觸時,實(shí)際接觸面積只占名義接觸面積的一小部分。微觀上認(rèn)為接觸是一個面上的一系列微小凸起與另外一個面上的微小凸起的相互抵觸作用?？蓪?個接觸面等效為一個理想光滑平面與另一個分布有微凸起的面。該簡化等效模型的當(dāng)量彈性模量可表示為

(2)

式中:E1、E2分別表示兩接觸物體的彈性模量;ν1、ν2分別表示2個接觸物體的泊松比。

(3)

(4)

1.2 切、法向耦合剛度的求解

微接觸面積a′的統(tǒng)計學(xué)分布[9]為

(5)

接觸剛度由接觸面間一系列微小凸起的彈-塑性變形引起。圖2為一微凸起模型,模型假設(shè)所有相接觸的凸起在接觸區(qū)域?yàn)椴煌霃降那蛐蝃7-9]。圖2中,r為微接觸面真實(shí)接觸半徑,可由a=πr2計算得出,a為實(shí)際接觸面積,即球形微凸起經(jīng)接觸壓縮后的接觸面積。對于一個圓形彈性微接觸,a′=2a[8-9],所以有a′=2πr2。

圖2 接觸微觀等效模型

若不計滑移產(chǎn)生的結(jié)構(gòu)拓?fù)渥冃?在接觸面有相對滑移的情況下,仍認(rèn)為接觸區(qū)域?yàn)閳A形。壓縮截面高度δ由W-M函數(shù)確定[8,10]

(6)

粗糙面上球形微凸起的半徑R可表示為R2=(R-δ)2+(r′)2,由于粗糙面的曲率半徑要遠(yuǎn)大于壓縮截面的高度δ,因此有(r′)2=2Rδ。根據(jù)赫茲接觸理論,微凸起的法向彈性變形力

(7)

單一微接觸對的法向接觸剛度kn可寫為

(8)

將式(6)、式(7)代入式(8),法向接觸剛度可寫為

(9)

根據(jù)Mindlin理論,切向剛度可寫為

kT=4G*r/(2-ν*)

(10)

式中:G*為材料的剪切模量;ν*為泊松比。

由以上各式,整體模型的接觸剛度[9]可寫為

(11)

整體模型的切向接觸剛度可寫為

(12)

由式(11)、式(12),每一個接觸點(diǎn)對的接觸剛度可用如下公式計算

(13)

式中:kni、kti分別表示第i個接觸點(diǎn)對的法向接觸剛度和切向接觸剛度;ai為第i個接觸點(diǎn)對的接觸面積。圖3所示為不同粗糙度下的切、法向接觸剛度值。

(a)法向接觸剛度

(b)切向接觸剛度圖3 接觸剛度與法向力的關(guān)系

2 求解扇區(qū)結(jié)構(gòu)振動的復(fù)約束法[4-5]

循環(huán)對稱結(jié)構(gòu)的單扇區(qū)有限元模型節(jié)點(diǎn)的位移可表示為

(14)

(15)

式中:n為節(jié)徑數(shù);N為組成循環(huán)對稱結(jié)構(gòu)的扇區(qū)數(shù)。需要注意的是,在結(jié)構(gòu)扇區(qū)有限元模型中,位于扇區(qū)一側(cè)的節(jié)點(diǎn)屬于該扇區(qū),而位于另一側(cè)的節(jié)點(diǎn)屬于相鄰扇區(qū),所以通常將有限元扇區(qū)模型稱為擴(kuò)充扇區(qū)。依據(jù)波動理論,扇區(qū)一側(cè)的節(jié)點(diǎn)位移可由另一側(cè)的節(jié)點(diǎn)位移表示,即

(16)

式(16)即為復(fù)約束條件,其中i為復(fù)數(shù)中的虛數(shù)單位,i2=-1。將式(16)代入式(14),得

(17)

(18)

(19)

3 耦合結(jié)構(gòu)振動特性求解方法

3.1 扇區(qū)有限元剛度矩陣及質(zhì)量矩陣的獲取

可通過修改Abaqus的inp文件進(jìn)行結(jié)構(gòu)扇區(qū)剛度矩陣及質(zhì)量矩陣的導(dǎo)出操作。用文本工具打開扇區(qū)模型的inp文件,在最后添加如下代碼命令:

*step

*matrix generate, stiffness, mass

*matrix output, stiffness, mass

*end step

代碼第2行指令用于生成剛度及質(zhì)量矩陣,第3行指令將這2個矩陣導(dǎo)出。返回Abaqus建立新作業(yè),在source項(xiàng)中選擇上述修改后的inp文件,提交后Abaqus將在工作目錄生成2個.mtx文件,分別代表剛度陣和質(zhì)量陣[6]。該.mtx文件包含5個列向量,第1、3列表示節(jié)點(diǎn)號,第2、4列為第1、3列的節(jié)點(diǎn)號對應(yīng)的自由度,第5列表示其對應(yīng)的剛度值。應(yīng)用MATLAB編寫代碼,將導(dǎo)出文件轉(zhuǎn)化為方陣或其稀疏格式。整理后的剛度陣是一個三角陣,這是因?yàn)閯偠染仃囀菍ΨQ的。質(zhì)量矩陣是維數(shù)為自由度3倍的對角矩陣。

3.2 節(jié)點(diǎn)尋找和歸類

需要找到3類節(jié)點(diǎn):①扇區(qū)主面節(jié)點(diǎn)和與其對應(yīng)的從面節(jié)點(diǎn);②結(jié)構(gòu)件與附加件結(jié)合部位的對應(yīng)耦合節(jié)點(diǎn);③求解時被約束的節(jié)點(diǎn)。Abaqus的inp文件中包含節(jié)點(diǎn)位置信息,使節(jié)點(diǎn)尋找和歸類變?yōu)榭赡堋Ｏ铝袀未a表示對節(jié)點(diǎn)的尋找和歸類:

1 導(dǎo)入.inp文件,將節(jié)點(diǎn)坐標(biāo)系轉(zhuǎn)化至柱坐標(biāo),形成矩陣C;

2 依據(jù)周向、徑向、軸向順序?qū)?jié)點(diǎn)進(jìn)行排列;

3 jj=1;kk=1;ll=1;mm=1;

4 for ii=節(jié)點(diǎn)1:最大節(jié)點(diǎn)

5 if節(jié)點(diǎn)周向位置-最小周向值<=tol

%tol為誤差限

6 left(jj,1)=C(ii,1);jj=jj+1;

%left為扇區(qū)主面節(jié)點(diǎn)集

7 if節(jié)點(diǎn)周向位置-最大周向值<=tol

8 right(kk,1)=C(ii,1);kk=kk+1;

%right為從面節(jié)點(diǎn)集

9 if節(jié)點(diǎn)徑向位置-軸孔半徑<=tol

10 n_bot(ll,1)=C(ii,1);ll=ll+1;

%n_bot為齒輪內(nèi)徑節(jié)點(diǎn)集

11 if節(jié)點(diǎn)軸向位置-耦合點(diǎn)軸向位置<=tol &&節(jié)點(diǎn)徑向位置-耦合點(diǎn)徑向位置<=tol

12

cont(mm,1)=C(ii,1);mm=mm+1;

%cont為耦合節(jié)點(diǎn)集

13end

在上面的偽代碼中:矩陣C為inp文件中節(jié)點(diǎn)信息轉(zhuǎn)換至柱坐標(biāo)系后的節(jié)點(diǎn)信息表;%后為注釋部分;左端數(shù)字為指令序號;第3條指令為各類節(jié)點(diǎn)的數(shù)目標(biāo)識;第4～12條為節(jié)點(diǎn)的尋找和歸類,其中第5～8條為扇區(qū)主、從界面節(jié)點(diǎn)尋找及存儲,第9、10條為約束節(jié)點(diǎn)尋找及存儲,第11、12條為耦合節(jié)點(diǎn)的尋找及存儲。

3.3 矩陣組裝及耦合剛度模擬

將組裝后的耦合結(jié)構(gòu)矩陣定義為2×2塊矩陣,其中將結(jié)構(gòu)扇區(qū)剛度矩陣記為K11,附加件扇區(qū)剛度矩陣記為K22,K12和K21為兩者間的耦合。用空間彈簧元模擬的耦合剛度,本質(zhì)上是在2個對應(yīng)節(jié)點(diǎn)相應(yīng)方向上的非零剛度系數(shù)。圖4為裝配了彈簧元的扇區(qū)擴(kuò)充剛度陣圖,2個扇區(qū)只在耦合節(jié)點(diǎn)處存在剛度耦合,其他節(jié)點(diǎn)處無耦合,應(yīng)補(bǔ)0。耦合節(jié)點(diǎn)對通過節(jié)點(diǎn)位置信息獲得(見3.2節(jié)),按式(13)賦值。

圖4 扇區(qū)擴(kuò)充剛度矩陣和彈簧元組裝

由于質(zhì)量間沒有耦合關(guān)系,故只需對質(zhì)量矩陣進(jìn)行擴(kuò)維,裝入附加件扇區(qū)質(zhì)量即可。

3.4 在主、從對稱面上添加復(fù)約束

循環(huán)對稱結(jié)構(gòu)的剛度矩陣和質(zhì)量矩陣在每個扇區(qū)建立的各自獨(dú)立的坐標(biāo)系下才具備塊循環(huán)特性,而有限元軟件是在整體直角坐標(biāo)系下導(dǎo)出剛度陣和質(zhì)量陣的。擴(kuò)充扇區(qū)的從面節(jié)點(diǎn)屬于下個扇區(qū),所以要對與從面節(jié)點(diǎn)相關(guān)的剛度陣做變換。擴(kuò)充扇區(qū)主面上的任意結(jié)點(diǎn)i的位移zi與從面上對應(yīng)的結(jié)點(diǎn)j的位移zj之間應(yīng)滿足如下循環(huán)對稱條件

(20)

式中:To為結(jié)點(diǎn)循環(huán)對稱變換矩陣;μ為每個輪齒所對應(yīng)的圓心角,μ=2π/N。上式是以z軸為軸向的,若結(jié)構(gòu)的軸向不同,變換矩陣也需要相應(yīng)調(diào)整。子結(jié)構(gòu)1的一般有限元特征值方程為

K1Z1=λM1Z1

(21)

由式(16)、式(17)和式(20)、式(21)可知,要建立基本扇區(qū)的復(fù)約束特征值方程(18),只需建立擴(kuò)充扇形子結(jié)構(gòu)的有限元模型,引入主、從面節(jié)點(diǎn)循環(huán)對稱條件,組裝擴(kuò)充扇形子結(jié)構(gòu)的剛度陣和質(zhì)量陣。復(fù)約束剛度陣和質(zhì)量陣[4-5]可寫為

(22)

4 求解實(shí)例

以弧齒錐齒輪-阻尼碗耦合結(jié)構(gòu)為例,根據(jù)上述方法,應(yīng)用MATLAB程序計算其振動特性。圖5為該耦合結(jié)構(gòu)的裝配關(guān)系:阻尼碗扣于錐齒輪背錐上,調(diào)節(jié)圓螺母和墊片為背錐提供不同的接觸力。齒輪扇區(qū)有限元模型(包含墊片扇區(qū))共有919個節(jié)點(diǎn),652個塊單元,如圖6所示。阻尼碗扇區(qū)有限元模型含189個節(jié)點(diǎn),80個塊單元,耦合后的扇區(qū)模型包含18個彈簧元。阻尼碗用錳鋼制作,厚度為1.5 mm,彈性模量E=210 GPa,泊松比ν=0.27。表1中的第3列為計算結(jié)果。

圖5 弧齒錐齒輪-阻尼碗裝配示意圖

圖6 齒輪扇區(qū)有限元模型

整體齒輪模型含59齒,54 221個節(jié)點(diǎn),38 468個塊單元;阻尼碗有7 434個節(jié)點(diǎn),4 720個塊單元;耦合模型含1 062個彈簧元。用Python編寫命令流將彈簧元添加至耦合位置,用Abaqus求解整體模型,結(jié)果列于表1的第4列。由于振型的正交性,零節(jié)徑之外的振型應(yīng)成對出現(xiàn)。由于網(wǎng)格及約束處理方式不同,在整體模型求解出的節(jié)徑振型中,同一振型對應(yīng)的共振頻率有微小偏差(16 Hz以內(nèi)),計算誤差時取其平均值。由表1可見,特征頻率誤差基本控制在1%內(nèi),說明本文提出的方法可以較好地模擬循環(huán)對稱耦合結(jié)構(gòu)的振動特性。

表1 2種模型計算的特征頻率對照表

Abaqus求解整體模型耗時152 s,應(yīng)用本文建立的扇區(qū)耦合模型求解程序(MATLAB編寫)耗時49 s,節(jié)省了約67%的求解時間,如果模型的扇區(qū)數(shù)量更多,則獲益將更大。

在表2的第2列中,振型代號短橫線前的數(shù)字表示節(jié)徑數(shù),后面的數(shù)字表示對應(yīng)的振型階數(shù);“碗”字表示該頻率對應(yīng)的振型以阻尼碗振動為主。圖7為對應(yīng)于表2中第12、13階頻率的振型,呈現(xiàn)出以阻尼碗振動為主的二節(jié)徑振動。

圖7 以阻尼碗振動為主的二節(jié)徑振動

5 結(jié) 論

(1)應(yīng)用W-M分形模型結(jié)合赫茲接觸理論和Mindlin理論,確定了耦合界面接觸剛度;根據(jù)結(jié)構(gòu)耦合后仍然呈現(xiàn)循環(huán)對稱這一特點(diǎn),用MATLAB開發(fā)出了利用結(jié)構(gòu)單扇區(qū)求解耦合結(jié)構(gòu)振動特性的計算流程。將該流程應(yīng)用于弧齒錐齒輪-阻尼碗耦合結(jié)構(gòu),并將所得的結(jié)果與耦合結(jié)構(gòu)的整體有限元模型計算結(jié)果進(jìn)行了對比,發(fā)現(xiàn)計算誤差基本都控制在1%之內(nèi),說明本文算法可以較好地模擬循環(huán)對稱耦合結(jié)構(gòu)的振動特性。同時,本文算法可以有效地節(jié)省存儲資源和計算時間,算例的求解時間由Abaqus的152 s縮短為49 s,節(jié)省了約67%。

(2)與阻尼環(huán)、阻尼襯套等較小附加件[2]不同,阻尼碗面積較大,且改變了齒輪的空間拓?fù)?因而振動分析結(jié)果增加了以阻尼碗振動為主的振型,在后續(xù)結(jié)構(gòu)減振計算中應(yīng)該注意到這一點(diǎn)。

[1] 晏礪堂, 朱梓根, 李其漢, 等. 高速旋轉(zhuǎn)機(jī)械振動 [M]. 北京: 國防工業(yè)出版社, 1994: 105-106.

[2] 曾亮, 李琳. 具有接觸結(jié)合面的篦齒封嚴(yán)組件振動特性分析 [J]. 航空動力學(xué)報, 2006, 21(5): 855-861. ZENG Liang, LI Lin. Vibration performance of contact system of labyrinth assembly [J]. Journal of Aerospace Power, 2006, 21(5): 855-861.

[3] THOMAS D L. Dynamics of rotationally periodic structures [J]. Int J Num Meth Engrg, 1979, 14(1): 81-102.

[4] RAMAMURTI V, BALASUBRAMANIAN P. Static analysis of circumferentially periodic structures with Potter’s scheme [J]. Computers & Structures, 1986, 22(3): 427-431.

[5] BALASUBRAMANIAN P, RAMAMURTI V. Free-vibration analysis of cyclic symmetric structures [J]. Comm in Appl Num Meth, 1991, 7(1): 131-139.

[6] Hibbitt, Karlsson and Sorensen, Inc. Abaqus 6.12 analysis user’s manual 2012: vol. 2 [M]. Pawtucket, USA: Abaqus Inc. Dassault Systems, 2012: 939-945.

[7] MAJUMDAR A, TIEN C L. Fractal characterization and simulation of rough surfaces [J]. Wear, 1990, 136(2): 313-327.

[8] YAN W, KOMVOPOULOS K. Contact analysis of elastic-plastic fractal surfaces [J]. J Appl Phys, 1998, 84(7): 3617-3624.

[9] JIANG S, ZHENG Y, ZHU H. A contact stiffness model of machined plane joint based on fractal theory [J]. J Tribol, 2010, 132: 011401.

[10]LIU Yalin, SHANGGUAN Bo, XU Zili. A friction contact stiffness model of fractal geometry in forced response analysis of a shrouded blade [J]. Nonlinear Dynamics, 2012, 70(3): 2247-2257.

(編輯 葛趙青)

Study on the Eigenvalue Problems of a Coupled Structure with Cyclic Symmetric Features

LIN Yanhu,ZHAO Ning,PENG Yanjun

(School of Mechanical Engineering, Northwestern Polytechnical University, Xi’an 710072, China)

Considerable amount of time and computing resources will be saved when the cyclic symmetric features are considered in the solution of eigenvalue problems if the components still have cyclic symmetric features after they were coupled. In this study, tangential and normal contact stiffness was solved by combining the Weierstrass-Mandelbrot (W-M) fractal geometry model and Herz contact theory. Then, the stiffness and mass matrices of each component were obtained by commercial finite element software, the total stiffness and mass matrices of the coupled structure were assembled by enlarging the obtained matrices, and the spring elements were used to simulate the contact stiffness of the coupled structure. The complex constraint condition was applied to the coupled structure, and the program for solving the eigenvalue problems by using only one sector model was developed with the software MATLAB. The program was used to analyze the eigenvector problem of a coupled aircraft spiral bevel gear-damper bowl structure, and the results were verified by the overall finite element model of the coupled structure. It is shown that the developed program could deal with the eigenvector problems with high efficiency and enough accuracy. The relative error was less than 1%, and 67% of computation time was saved, meanwhile, the memory resource was also saved. In addition, it is shown that the eigenvalues were obtained, accompanied with the eigenvectors revealing the vibration of the damper bowl. So sufficient attention should be paid in the following calculation including frequency response or vibration alleviation.

coupled structure; cyclic symmetric structure; vibration; eigenvalue problem; complex constraint

2015-11-01。 作者簡介:藺彥虎(1983—),男,博士生;趙寧(通信作者),男,教授。 基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目(51075328)。

時間:2016-04-28

10.7652/xjtuxb201607011

TH113

A

0253-987X(2016)07-0068-07

網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20160428.2222.004.html