999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

滾動軸承-轉子系統的分岔分析與多吸引子共存

2022-04-12 00:00:00安慧寧金花呂小紅王昕
機械強度 2022年3期

摘要以滾動軸承支撐下的不平衡轉子系統為研究對象,運用變步長 Runge-Kutta 法進行數值積分獲取轉子系統的動力學響應。根據增速和減速不同工況下的分岔圖、相圖和 Poincare 映射圖分析系統共存吸引子的分岔,揭示其吸引域隨系統參數變化的演變過程。結果表明:隨著轉速的變化,系統在發生 Hopf 分岔、跳躍分岔和倍化分岔時,會出現吸引子共存的現象。跳躍分岔會導致其吸引域的拓撲結構發生突變,而 Hopf 分岔和倍周期分岔對其吸引域的影響較小。該研究結果可為系統在不同轉速下運行時提供指導,為滾動軸承-轉子系統的平穩運行提供理論依據。

關鍵詞滾動軸承轉子系統跳躍分岔共存吸引子吸引域

中圖分類號 TH113

AbstractTaking the unbalanced rotor system supported by rolling bearings as the research object , the variable step sizeRunge-Kutta method isused fornumericalintegrationtoobtainthedynamicresponseof therotorsystem. Accordingtothebifurcation diagrams , phase portraits and Poincare maps under different conditions of increasing and decelerating , the bifurcation of coexisting attractor is analyzed , and the evolution process of its attraction region with the system parameters was revealed . The results show that when Hopf bifurcation , jump bifurcation and doubling bifurcation occur in the system with the change of rotating speed ,there will be the coexistence of attractors . The jump bifurcation will lead to a sudden change in the topological structure of the basin of attraction , while the Hopf bifurcation and the doubling bifurcation have little effect on the basin of attraction. The research results can provide guidance for the operation of the system at different speeds , and provide theoretical basis for the smooth operation of the rolling bearing rotor system.

Key wordsBall bearing;Rotor system;Jump bifurcation;Coexistence attractors;Basin of attraction

Corresponding author:JIN Hua , E-mail : jinh@ mail. lzjtu.com

The project supported by the National Natural Science Foundation of China ( No .12062008), and the Gansu Science and Technology Planning Project ( No .20YF8WA043,20JR5RA410).

Manuscript received 20201117 in revised form 20210511.

引言

滾動軸承支撐的轉子系統廣泛應用于火箭發動機、航空燃氣渦輪機及其他旋轉機械系統,其動力學特性的研究得到人們極大的關注[1]。近年來,對于滾動軸承系統的動力學特性研究已有很多成果。 Tiwari M 等[2]研究了不平衡力、軸承游隙、VC 振動和非線性 Hertz 接觸力對 Jeffcott 轉子非線性動力學特性的影響。 Lioulios A N 等[3]通過研究變轉速對系統動態特性的影響,發現轉子的速度波動是系統動態行為的控制參數。袁茹等[4]研究了系統響應隨轉速的變化趨勢,發現軸承的徑向游隙是決定系統動態響應的一個關鍵參數。陳果[5]2773-2778建立了含故障不平衡轉子-滾動軸承動力學模型,發現了通往混沌的倍周期分岔的途徑。張偉剛等[6]研究了具有負游隙的機床主軸-滾動軸承系統的非線性動態特性和穩定性。白長青等[7]研究了滾動軸承平衡轉子系統在不同軸承內的間隙量、不同轉速下系統的穩定性及其分岔特性和混沌。陳恩利等[8]建立了含滾動軸承支撐松動故障的動力學模型,研究了系統分岔和混沌運動等復雜動力學現象及變化規律。呂運等[9]建立了滾動軸承轉-子-基座的動力學模型,分析了不同參數下的動力學特性。馮吉路等[10]建立了考慮彈性主軸剛度、阻尼和滾動軸承非線性接觸力的十自由度數控機床主軸非線性振動分析模型。李默等[11]建立了含橫傾角的非對稱滾動軸承-轉子動力學模型,研究了橫傾角對系統的影響。

上述文獻建立了不同自由度的軸承轉子系統動力學模型,都是側重于系統的分岔行為及周期運動通向混沌的路徑分析,不能全面地揭示系統的動力學行為及全局特性,對于系統中共存周期解及其分岔演化過程鮮有提及。因此,本文使用文獻[5]2773-2778建立的具有質量偏心的轉子系統非線性動力學模型,分析滾動軸承轉子系統隨轉速變化的運動狀態;根據轉速增大和減小兩種工況下的分岔圖,研究系統周期運動的共存分岔模式,并結合吸引域、相圖和 Poincare 映射圖揭示系統中共存吸引子的吸引域隨轉速變化的演化規律。

1 滾動軸承支撐下的轉子動力學模型

1.1滾動軸承模型

本文所選取的滾動軸承如圖1所示,外滾道半徑為 Ro=63.9 mm、內滾道半徑為 Ri=40.1 mm、滾珠半徑 rb=5.59 mm、滾珠個數 Nb=8個。軸承外圈和內圈的旋轉角速度分別用ω o 、ω i 表示。滾動軸承外圈固定在剛性體上,所以滾珠與外圈接觸點的線速度 vo=0;內圈固定在旋轉軸上,故內圈線速度 vi =ω i Ri =ωRi ,則保持架的角速度ω c =。將各個滾子上的接觸力求和,同時考慮軸承阻尼,可得到前、后軸承沿 x 、y 方向總的接觸力

式中,θ i =ω c t +( i -1)表示第 i 個滾珠的角位置。

γ為滾動軸承的徑向內部間隙,γ=5 μm ,則第 i 個滾子上接觸點的彈性變形 r θi =xcosθ i +ysin θ i -γ,下標“+”表示方括號內的值為非負型。接觸點的接觸剛度系數kb 為

kb =(2)

式中,E 為彈性模量,1/M 為泊松比, p 為接觸單元的曲率之和。根據參考文獻[12],取 E=2.1×1011 N/m2,1/M=0.3, p =266 m-1,代入式(2)可得接觸剛度 kb=13.34×109 N/m2/3。

滾動軸承在旋轉時,由于滾動體依次通過載荷區,在徑向載荷作用下會產生有滾動體和無滾動體兩種情況,二者的接觸剛度不同,因此會引起接觸變形的差異,導致參數激振,稱為 VC(Varying Compliance)振動現象。其頻率可表示為滾動體公轉角速度與滾動體個數的乘積

fvc =ω c Nb(3)

1.2系統動力學模型及運動方程

本文考慮滾動軸承-轉子系統[5]2773-2778。為了方便分析,將該系統簡化為三個質量單元,并且只考慮質量單元在 x 和 y 方向的平動,其簡化后的六自由度系統如圖2所示。 O 1、O2、O3分別為軸承的幾何中心、轉子的幾何中心和轉子質心;mp 為轉子在圓盤處的等效集中質量, mL 和 mR 為轉子在左右兩端軸承處的轉子集中質量;cb 和 cp 分別為轉子在軸承處和圓盤處的阻尼系數;k 為彈性軸剛度;fxL 和 fyL 為左端軸承的支撐反力,fxR 和fyR為右端軸承的支撐反力。

根據牛頓第二定律,得到滾動軸承-轉子系統的動力學微分方程

式中,xp 、yp 分別為轉子圓盤處 x 方向和 y 方向的位移;xR 、yR 分別為轉子右軸承處 x 方向和 y 方向的位移;xL 、yL 分別為轉子左軸承處 x 方向和 y 方向的位移;e 為系統的質量偏心。

為使該類模型的分析具有一般性,取無量綱參數

得無量綱動力學微分方程

其中

2 周期運動的分岔

現有對轉子系統動力學的研究中,都是側重于系統的分岔行為及周期運動通向混沌的路徑分析。這種分析方法不能全面地揭示系統的動力學行為及全局特性,而且非線性系統的振動響應對系統的初始條件比較敏感,為了能夠對工程實際提供更多的理論依據,本文基于分岔參數遞增和遞減的不同工況,研究系統周期運動的共存分岔模式及共存吸引子的吸引域演化,揭示系統的分岔特點及全局動力學行為。

本文選取的滾動軸承-轉子系統基本參數如下:mR= mL =4.0 kg , mp =32.1 kg , e =0.01 mm , cb =1050 N .s/m ,cp=2100 N . s/m , k=10×107 N/m ,計算圓盤位移 xp 隨轉速ω增大的分岔圖如圖3a 所示。圖3b 為 xp 隨轉速ω減小和增大的合成分岔圖,其中深色為ω減小時的分岔圖,淺色為ω增大時的分岔圖。由圖3a 可見,當轉速ω∈[945,2037]區間內,系統表現為周期一運動。減小ω,周期一運動經倍周期分岔產生周期二運動,然后經 Hopf 分岔產生周期二運動的擬周期運動。周期二運動及其擬周期運動的 Poincare 映射圖分別如圖4a 和圖4b 所示。當轉速較低時,由于旋轉不平衡引起的振動響應比較弱,而軸承總體剛度周期變化引起的參數激勵較強,因此,當ω小于800 rad/s 時,系統會出現混沌運動和周期一運動的擬周期運動交替出現的現象。圖4c 和圖4d 分別為ω=700 rad/s 和ω=450 rad/s 時混沌運動和擬周期運動的 Poincare 映射圖。而當增大ω至ω= 2037 rad/s 時,周期一運動發生跳躍分岔產生周期三運動,然后在ω= 2257 rad/s 時又經跳躍分岔產生周期二運動。圖5a 和圖5b 分別為跳躍分岔前后周期一運動和周期三運動的相圖和 Poincare 映射圖;圖5c 和圖5d 分別為跳躍分岔前后周期三運動和周期二運動的相圖和 Poincare 映射圖。因為非線性振動系統的振幅響應是振動工程在實際應用中最為關心的問題之一,且當系統的振幅超過一定限度,往往會導致系統結構的破壞甚至發生安全事故。由圖5可見,在跳躍分岔前后,圓盤的振動幅值發生了突變。在周期一運動的區間內,系統主要受到轉子不平衡激勵所引起的強迫振動,其位移振幅遠大于周期三運動時的振幅。因此,當轉子系統處于加速工況時,為保持系統的平穩加速應該盡快通過跳躍分岔區。

由圖3b 可見,周期一運動隨轉速增大時的分岔過程是不可逆的,因此在ω∈[1750,2260]區間內出現了周期一 P 1 a 和周期一 Q 1 a 、周期三 P3 a 和周期一 Q 1 a 、周期三 P3 a 和周期二 Q2a的共存區。在ω∈[2700,2990]區間內,系統存在周期二 P2b 和周期二 Q2b ,P2b 和周期二的擬周期運動,以及 P2b 和混沌運動的共存區。且在該區間內,隨著ω的減小,系統發生周期二運動的逆 Hopf 分岔,其中圖6 所示的 Poincare 映射圖描述了混沌運動經擬周期運動的環面變形、擬周期運動轉遷為周期二運動的過程。如圖3b 所示,繼續減小ω,周期二 Q2b 發生跳躍分岔轉遷為 P2b ,且隨ω減小的 Hopf 分岔過程是不可逆的。當ω增大時,周期二 P2b 經跳躍分岔直接轉遷為混沌。

由圖3 a 可見,在ω∈[3234,3513]區間內,系統表現為周期一運動。當ω減小時,周期一運動的分岔過程為:周期一運動→周期二運動→擬周期運動→環面變形→混沌,見圖7。當ω增大時,周期一運動經 Hopf 分岔產生擬周期運動。

3 吸引域的演化

由第2節的分析可知,隨著轉速的變化,系統發生 Hopf 分岔和跳躍分岔等動力學行為。由于跳躍分岔的不可逆性,系統在跳躍分岔點鄰域內存在共存分岔模式和吸引子共存的現象。當多吸引子共存時,系統的初始條件對動力學響應的影響很大,因此本節在系統相空間中選取一個考察區域 H={-20.0≤ Xp ≤20.0,-20.0≤ X(.)≤20.0},通過計算共存吸引子在 H 內的吸引域揭示系統的全局動力學行為及吸引域的演化過程。由圖3b 可知,在ω∈[1750,2260]內,系統中存在兩個不同的周期一吸引子 P1 a 和 Q1 a 共存、周期三 P3 a 和周期一 Q1 a 共存及周期三 P3 a 和周期二 Q2a 共存的現象,共存吸引子的吸引域如圖8所示。當ωlt;1750 rad/s 時,系統只存在穩定的周期一吸引子 P1 a 。當ω增加至ω= 1750 rad/s 時,P1 a 在部分初值范圍內發生周期跳躍,在考察區域 H 中產生新的周期一吸引子 Q1 a 。因此當ωgt;1750 rad/s時,周期一吸引子 P1 a 和 Q1 a 共存。圖8a 給出了ω= 1850 rad/s 時共存周期一吸引子 P1 a 和 Q1 a 的吸引域,圖中周期一吸引子 P1 a 未在考察區域 H 內,因此未標注,淺色區域為其吸引域;“×”號代表周期一吸引子 Q1 a ,深色區域為其吸引域。由圖可見, P1 a 的淺色吸引域中心地帶被 Q1 a 的深色吸引域侵蝕,即 P1 a 在部分初值范圍內發生周期跳躍轉遷為 Q1 a ,但是淺色區域的面積遠大于深色區域的面積,因此 P1 a 的穩定性最強。當ω= 1950 rad/s 時,P1 a 和 Q1 a 的吸引域如圖8b 所示,對比圖8a ,可見隨著轉速ω的增大, P1 a 的吸引域被 Q1 a 的吸引域大量侵蝕, Q1 a 的吸引域面積逐漸擴大,其穩定性逐漸增強。圖9a 為 P1 a 和 Q1 a 的相圖和 Poincare 映射圖,由圖可見, P1 a的振動幅值遠大于 Q1 a 的振動幅值,故 P1 a 對系統振動極為不利,容易引起系統結構的破壞,在實際運動中應避免此類運動狀態的產生。

進一步增加ω,當ω=2037 rad/s時,周期一吸引子 P1 a 經跳躍分岔產生周期三吸引子 P3 a ,因此當ω∈[2037,2165]時,周期三吸引子 P3 a 和周期一吸引子 Q1 a 共存。圖8c 為ω=2060 rad/s的兩個共存吸引子的吸引域,“+”代表周期三吸引子 P3 a ,淺色區域為其吸引域;“×”代表周期一吸引子 Q1 a ,深色區域為其吸引域。由圖可見, P3 a 的淺色吸引域和 Q1 a 的深色區域相互嵌套并纏繞在一起,此時系統在該轉速下兩個吸引子的全局動力學特性都不穩定,只要初值稍有變化就會使系統運動到不同的周期軌道上。圖8d 為ω=2150 rad/s時 P3 a 和 Q1 a 的吸引域,對比圖8c 可見,隨著ω增加, Q1 a 的吸引域面積不斷增大,邊界趨于光滑,P3 a 的吸引域被 Q1 a 的吸引域不斷侵蝕,面積減小,吸引子 Q1 a 的穩定性增強,圖9b 為ω= 2150 rad/s 時共存吸引子 P3 a 和 Q1 a 的相圖和 Poincare 映射圖。

增大ω至ω= 2165 rad/s 時,周期一吸引子 Q1 a 經倍周期分岔產生周期二吸引子 Q2a ,并與周期三吸引子 P3 a 共存。圖8e 為ω= 2180 rad/s 時共存吸引子 P3 a 和 Q2a 的吸引域,“+”代表周期三吸引子 P3 a ,淺色區域為其吸引域;“×”代表周期二吸引子 Q2a ,深色區域為其吸引域。由圖可見, Q2a 和 P3 a 的吸引域纏繞程度低,分布相對集中,但是明顯深色吸引域的面積大于淺色吸引域的面積,故周期二吸引子 Q2a 的穩定性比周期三吸引子 P3 a 的穩定性強。結合圖8e 和圖8f 可知,繼續增大ω, Q2a 的深色吸引域面積不斷增加,P3 a 的吸引域被進一步侵蝕,兩共存吸引子的吸引域在 H 的右半側互相交替呈輻射狀分布,分界面趨于光滑,纏繞度繼續降低。圖9c 為ω= 2220 rad/s 時共存吸引子 P3 a 和 Q2a 的相圖和 Poincare 映射圖。

由圖8吸引域的演化過程可知,跳躍分岔對于系統全局動力學特性的影響較大,會引起吸引域拓撲結構的突變,而周期倍化分岔對共存吸引子的吸引域影響較小。

如圖3b 所示,當ω∈[2700,2990]范圍內時,隨著ω的減小,系統從混沌運動退化為擬周期運動;進一步減小ω,擬周期運動又退化為周期二運動 Q2b ;最后 Q2b 經跳躍分岔轉遷為周期二 P2b 。而ω在此范圍增大時,系統只存在周期二吸引子 P2b 。當ω =2805 rad/s 時,系統中兩個不同的周期二吸引子 P2b 和 Q2b 共存,其吸引域如圖10a 所示。“+”代表周期二吸引子 P2b ,淺色區域為其吸引域;“×”代表周期二吸引子 Q2b ,深色區域為其吸引域。 P2b 和 Q2b 的吸引域相互嵌套并纏繞在一起,但是淺色吸引域的面積相對集中,其對應的吸引子 P2b 穩定性較強。圖10b 為ω= 2850 rad/s時 P2b 和 Q2b 的吸引域,與圖10a 相比,隨著ω的增加, Q2b 的吸引域面積逐漸增大,在考察區域 H 內的左半部分其吸引域的分布也變得相對集中,P2b 被 Q2b 進一步侵蝕,此時 Q2b 的穩定性強于 P2b 的穩定性。圖11a 為共存吸引子 P2b 和 Q2b的相圖和 Poincare 映射圖。

繼續增大ω,由圖3b 可見,周期二吸引子 P2b 和周期二的擬周期吸引子共存,圖10c 和圖11b 分別為ω=2900 rad/s 時 P2b 和擬周期吸引子共存的吸引域圖、相圖和 Poincare 映射圖,圖中“+”代表周期二吸引子 P2b ,淺色區域為其吸引域;深色區域為擬周期吸引子的吸引域。圖10d 和圖11c 分別為ω=2980 rad/s時 P2b 和混沌吸引子共存的吸引域、相圖和 Poincare 映射圖,“+”代表周期二吸引子 P2b ,淺色區域為其吸引域;深色區域為混沌吸引子的吸引域。對比圖10c 和圖10d 可知,隨著ω的增大,P2b 的淺色吸引域從左向右不斷被侵蝕,面積逐漸減小。深色區域的面積逐漸集中,擬周期運動逐漸轉遷為混沌運動,混沌吸引子的穩定性強于周期二吸引子 P2b 。進一步增大ω,周期二吸引子 P2b 經跳躍分岔進入混沌。

由以上分析可知,當系統中共存吸引子的吸引域相互嵌套纏繞在一起時,系統對初值條件極為敏感,使得系統的穩定性降低;隨著轉速ω的增大,各吸引子的吸引域不斷集中,且由左側向右側逐漸擴散,使得系統各吸引子的穩定性有所增強。

4 結論

本文考慮滾動軸承 Hertz 接觸力和具有轉子偏心的動力學模型,研究了系統在轉速增大和減小兩種情況下周期運動的共存分岔模式,分析了共存吸引子吸引域的演化過程。得出以下結論:

1)系統存在 Hopf 分岔、跳躍分岔和倍化分岔等多種形式。由于變柔度 VC 振動和不平衡力導致的多頻激勵產生了多種形式的 Hopf 分岔。隨著轉速的變化,Hopf 分岔會產生周期二運動的擬周期運動;或發生逆 Hopf 分岔,使混沌運動經鎖相、擬周期運動轉遷為周期二運動。

2)跳躍分岔可導致周期一運動轉遷為周期三運動,或者直接由周期運動轉遷為混沌運動。且跳躍分岔前后,轉子的振動幅值突變,在臨界轉速附近,周期一運動振幅遠大于周期三運動時的振幅,對系統振動極為不利,容易引起系統結構的破壞,在實際運動中應合理控制系統轉速快而平穩的通過該區域。

3)系統在分岔參數遞增和遞減的不同工況下存在多吸引子共存的現象。隨著轉速的增大,新出現的吸引子對共存吸引子的侵蝕會逐漸增強,穩定性也逐漸增強。當共存吸引子的吸引域互相纏繞嵌套時,系統對初值較為敏感,各吸引子的穩定性較低,全局動力學特性復雜。當系統發生跳躍分岔時,對系統的全局動力學特性影響較大,會導致考察區域 H 內吸引域的拓撲結構發生突變。而倍化分岔和 Hopf 分岔對吸引域的拓撲結構影響較小。

參考文獻(References)

[1] 張耀強,陳建軍,唐六丁,等.滾動軸承-轉子系統非線性參數、強迫聯合振動[ J].機械強度,2009,31(6):871-875.

ZHANGYaoQiang , CHENJianJun , TANGLiuDing , etal . Nonlinearvibrationsofarollingbearing-rotorsystemsubjectto parametrical andexternalexcitations [ J ]. JournalofMechanical Strength ,2009,31(6):871-875( In Chinese).

[2]Tiwari M , Gupta K , Prakash O . Dynamic response of an unbalancedrotorsupportedonballbearings [ J ]. JournalofSoundand Vibration ,2000,238(5):757-779.

[3]Lioulios A N , Antoniadis I A . Effect of rotational speed fluctuationsonthedynamicbehaviourof rollingelementbearingswithradial clearances [ J]. International Journal of Mechanical Sciences ,2006(48):809-829.

[4] 袁茹,趙凌燕,王三民.滾動軸承轉子系統的非線性動力學特性分析[ J].機械科學與技術,2004(10):1175-1177.

YUANRu , ZHAOLinYan , WANGSanMing. Analysisofthe nonlinear dynamic behaviors of a rolling bearing-rotor system [ J]. MechanicalScienceandTechnologyforAerospaceEngineering ,2004(10):1175-1177( In Chinese).

[5] 陳果.滾動軸承支承下的不平衡轉子系統非線性動力響應分析[ J].中國機械工程,2007(23):2773-2778.

CHEN Guo . Nonlinear dynamic response analysis of an unbalanced rotor supported on ball bearing [ J]. China Mechanical Engineering ,2007(23):2773-2778( In Chinese).

[6] 張偉剛,高尚晗,龍新華,等.機床主軸-滾動軸承系統非線性動力學分析[ J].振動與沖擊,2008,27(9):72-75.

ZHANGWeiGang , GAOShangHan , LONGXinHua , etal . Nonlinear analysis for a machine-tool spindle system supported with ball bearings [ J]. Journal of Vibration and Shock ,2008(9):72-75( In Chinese).

[7] 白長青,許慶余,張小龍.考慮徑向內間隙的滾動軸承平衡轉子系統的非線性動力穩定性[ J ].應用數學和力學,2006(2):159-170.

BAI ChangQing , XU QingYu , ZHANG XiaoLong. Nonlinear stability of balanced rotor due to the effect of ball bearing internal clearance[ J]. Applied Mathematics and Mechanics ,2006(2):159-170( In Chinese).

[8] 陳恩利,王洪禮,何田.轉子系統支承松動的分岔特性與故障診斷[ J].機械強度,2007(3):387-389.

CHENEnLi , WANGHongLi , HETian . Bifurcationandfault diagnosis ofrotorsystemwithsupportloosening [ J ]. Journalof Mechanical Strength ,2007(3):387-389( In Chinese).

[9] 呂運,童大鵬,田野,等.滾動軸承-轉子系統動力學建模與仿真分析[ J].機械強度,2015,37(6):1178-1185.

Lü Yun , TONG DaPeng , TIAN Ye , et al . Dynamic modeling and simulation analysis of a rolling bearing-rotor system [ J]. Journal of Mechanical Strength ,2015,37(6):1178-1185( In Chinese).

[10] 馮吉路,汪文津,田越.數控機床主軸系統動力學建模與頻率特征分析[ J].組合機床與自動化加工技術,2017(8):29-32.

FENG JiLu , WANG WenJin , TIANYue . Dynamicsmodeling and frequency characteristic analysis for a machine spindle system [ J]. Modular Machine Tool Automatic Manufacturing Technique ,2017(8):29-32( In Chinese).

[11] 李默,劉永葆,王強.橫傾角影響下非對稱支撐轉子的非線性動力學研究[ J].機械強度,2020,42(2):293-298.

LI Mo , LIU YongBao , WANG Qiang. Study on nonlinear dynamics of asymmetrical supported rotor under influence of heeling angle [ J]. JournalofMechanicalStrength ,2020,42(2):293-298 ( In Chinese).

[12]Harsha SP , Kankar PK 1. Stabilityanalysisof arotorbearingsystemduetosurfacewavinessandnumberofballs [ J ]. International Journal of Mechanical Sciences ,2004,46(7):1057-1081.

主站蜘蛛池模板: 这里只有精品免费视频| 日韩无码黄色网站| 高清无码一本到东京热| 新SSS无码手机在线观看| 国产超薄肉色丝袜网站| 在线精品亚洲一区二区古装| 亚洲精选高清无码| 网久久综合| 久久国产香蕉| 视频一本大道香蕉久在线播放| 波多野一区| 成人亚洲视频| 亚洲综合色吧| 91年精品国产福利线观看久久| 中文无码影院| 99视频在线免费观看| 久久免费观看视频| 伊人欧美在线| 日韩精品一区二区三区免费在线观看| 99成人在线观看| 在线综合亚洲欧美网站| 91美女视频在线| 草草影院国产第一页| 日韩欧美综合在线制服| 亚洲一道AV无码午夜福利| 国产特级毛片| 成人福利一区二区视频在线| 国产人妖视频一区在线观看| 国产精品亚洲αv天堂无码| 激情综合五月网| 日本免费一级视频| www.亚洲色图.com| 久久天天躁夜夜躁狠狠| 日本爱爱精品一区二区| 亚洲国产天堂久久综合226114| 精品成人免费自拍视频| 国产一区亚洲一区| 国产在线自揄拍揄视频网站| 免费Aⅴ片在线观看蜜芽Tⅴ| 亚洲精品视频在线观看视频| 青青国产视频| a亚洲天堂| 色综合久久久久8天国| 亚洲av无码久久无遮挡| 少妇精品网站| 亚洲中文字幕手机在线第一页| 色偷偷一区| 久久午夜夜伦鲁鲁片无码免费| 国产制服丝袜91在线| 97se亚洲综合不卡| 亚洲三级网站| 蜜桃视频一区二区三区| 在线观看91精品国产剧情免费| 色呦呦手机在线精品| 青青草原国产| 欧美亚洲国产一区| 中文字幕日韩欧美| 天堂岛国av无码免费无禁网站| 中文字幕永久在线观看| 免费观看三级毛片| 国内精品手机在线观看视频| 人妖无码第一页| 国产欧美中文字幕| 午夜限制老子影院888| 毛片免费在线视频| 免费看美女自慰的网站| 国产91小视频在线观看| 亚洲一区二区三区国产精品 | 国产成人一级| 亚洲国产在一区二区三区| 特级aaaaaaaaa毛片免费视频 | 超薄丝袜足j国产在线视频| 久久成人国产精品免费软件 | 国产正在播放| 国产特一级毛片| 尤物亚洲最大AV无码网站| 福利在线一区| 欧美色亚洲| 粗大猛烈进出高潮视频无码| 免费无码又爽又黄又刺激网站| 毛片基地美国正在播放亚洲| 免费AV在线播放观看18禁强制|