具有循環對稱特點的耦合結構振動特性分析

藺彥虎,趙寧,彭艷軍

(西北工業大學機電學院,710072,西安)

?

具有循環對稱特點的耦合結構振動特性分析

藺彥虎,趙寧,彭艷軍

(西北工業大學機電學院,710072,西安)

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

耦合結構;循環對稱結構;振動;特征值;復約束

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

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

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

1 耦合剛度計算

1.1 W-M分形模型

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

(1)

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

圖1 不同分形維數的分形表面輪廓

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

(2)

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

(3)

(4)

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

微接觸面積a′的統計學分布[9]為

(5)

接觸剛度由接觸面間一系列微小凸起的彈-塑性變形引起。圖2為一微凸起模型,模型假設所有相接觸的凸起在接觸區域為不同半徑的球形[7-9]。圖2中,r為微接觸面真實接觸半徑,可由a=πr2計算得出,a為實際接觸面積,即球形微凸起經接觸壓縮后的接觸面積。對于一個圓形彈性微接觸,a′=2a[8-9],所以有a′=2πr2。

圖2 接觸微觀等效模型

若不計滑移產生的結構拓撲變形,在接觸面有相對滑移的情況下,仍認為接觸區域為圓形。壓縮截面高度δ由W-M函數確定[8,10]

(6)

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

(7)

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

(8)

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

(9)

根據Mindlin理論,切向剛度可寫為

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

(10)

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

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

(11)

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

(12)

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

(13)

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

(a)法向接觸剛度

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

2 求解扇區結構振動的復約束法[4-5]

循環對稱結構的單扇區有限元模型節點的位移可表示為

(14)

(15)

式中:n為節徑數;N為組成循環對稱結構的扇區數。需要注意的是,在結構扇區有限元模型中,位于扇區一側的節點屬于該扇區,而位于另一側的節點屬于相鄰扇區,所以通常將有限元扇區模型稱為擴充扇區。依據波動理論,扇區一側的節點位移可由另一側的節點位移表示,即

(16)

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

(17)

(18)

(19)

3 耦合結構振動特性求解方法

3.1 扇區有限元剛度矩陣及質量矩陣的獲取

可通過修改Abaqus的inp文件進行結構扇區剛度矩陣及質量矩陣的導出操作。用文本工具打開扇區模型的inp文件,在最后添加如下代碼命令:

*step

*matrix generate, stiffness, mass

*matrix output, stiffness, mass

*end step

代碼第2行指令用于生成剛度及質量矩陣,第3行指令將這2個矩陣導出。返回Abaqus建立新作業,在source項中選擇上述修改后的inp文件,提交后Abaqus將在工作目錄生成2個.mtx文件,分別代表剛度陣和質量陣[6]。該.mtx文件包含5個列向量,第1、3列表示節點號,第2、4列為第1、3列的節點號對應的自由度,第5列表示其對應的剛度值。應用MATLAB編寫代碼,將導出文件轉化為方陣或其稀疏格式。整理后的剛度陣是一個三角陣,這是因為剛度矩陣是對稱的。質量矩陣是維數為自由度3倍的對角矩陣。

3.2 節點尋找和歸類

需要找到3類節點:①扇區主面節點和與其對應的從面節點;②結構件與附加件結合部位的對應耦合節點;③求解時被約束的節點。Abaqus的inp文件中包含節點位置信息,使節點尋找和歸類變為可能。下列偽代碼表示對節點的尋找和歸類:

1 導入.inp文件,將節點坐標系轉化至柱坐標,形成矩陣C;

2 依據周向、徑向、軸向順序對節點進行排列;

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

4 for ii=節點1:最大節點

5 if節點周向位置-最小周向值<=tol

%tol為誤差限

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

%left為扇區主面節點集

7 if節點周向位置-最大周向值<=tol

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

%right為從面節點集

9 if節點徑向位置-軸孔半徑<=tol

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

%n_bot為齒輪內徑節點集

11 if節點軸向位置-耦合點軸向位置<=tol &&節點徑向位置-耦合點徑向位置<=tol

12

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

%cont為耦合節點集

13end

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

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

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

圖4 扇區擴充剛度矩陣和彈簧元組裝

由于質量間沒有耦合關系,故只需對質量矩陣進行擴維,裝入附加件扇區質量即可。

3.4 在主、從對稱面上添加復約束

循環對稱結構的剛度矩陣和質量矩陣在每個扇區建立的各自獨立的坐標系下才具備塊循環特性,而有限元軟件是在整體直角坐標系下導出剛度陣和質量陣的。擴充扇區的從面節點屬于下個扇區,所以要對與從面節點相關的剛度陣做變換。擴充扇區主面上的任意結點i的位移zi與從面上對應的結點j的位移zj之間應滿足如下循環對稱條件

(20)

式中:To為結點循環對稱變換矩陣;μ為每個輪齒所對應的圓心角,μ=2π/N。上式是以z軸為軸向的,若結構的軸向不同,變換矩陣也需要相應調整。子結構1的一般有限元特征值方程為

K1Z1=λM1Z1

(21)

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

(22)

4 求解實例

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

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

圖6 齒輪扇區有限元模型

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

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

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

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

圖7 以阻尼碗振動為主的二節徑振動

5 結 論

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

(2)與阻尼環、阻尼襯套等較小附加件[2]不同,阻尼碗面積較大,且改變了齒輪的空間拓撲,因而振動分析結果增加了以阻尼碗振動為主的振型,在后續結構減振計算中應該注意到這一點。

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

[2] 曾亮, 李琳. 具有接觸結合面的篦齒封嚴組件振動特性分析 [J]. 航空動力學報, 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—),男,博士生;趙寧(通信作者),男,教授。 基金項目:國家自然科學基金資助項目(51075328)。

時間:2016-04-28

10.7652/xjtuxb201607011

TH113

A

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

網絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20160428.2222.004.html