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

鼠籠彈性支承干摩擦減振試驗研究和參數識別

2023-12-08 02:26:28宋譚陳九如李浦袁奇
西安交通大學學報 2023年11期
關鍵詞:模型

宋譚,陳九如,李浦,袁奇

(西安交通大學能源與動力工程學院,710049,西安)

隨著大推重比、高轉速航空發動機的發展,柔性轉子在起動、加速、停車等過程中必須多次通過臨界轉速。鼠籠式彈性支承是典型的航空發動機支承結構,其剛度對轉子振動特性影響顯著[1]。若設計不當,可能引起發動機轉子振幅過大造成轉子斷裂,產生嚴重的后果。因此,開展鼠籠彈性支承剛度特性和減振研究具有重要意義[2]。

國內外學者針對鼠籠彈性支承結構的研究主要集中在剛度分析和試驗參數識別兩個方面。徐方程和劉占生開展了不同結構鼠籠支承剛度試驗研究,發現了鼠籠條數、籠條長度和籠條面積對其剛度的影響規律[3]。彭京徽等針對鼠籠彈性支承簡化模型,采用有限元方法開展了可靠性分析和降低應力集中方法的研究,表明倒圓角可以有效減少應力集中,而應力槽和圓孔對降低應力影響不大[4-6]。馮國全等基于參數化設計方法,考慮籠條非對稱彎曲特性的鼠籠式彈性支承柔度和應力的參數化方法,使用遺傳算法對彈性支承進行優化,大幅降低疲勞應力[7]。王洋洲等針對折返式鼠籠彈性支承剛度問題開展了理論分析和試驗驗證[8-9]。孫彥博等分析了網格層數對鼠籠彈性支承剛度計算的影響[10]。唐瑞等基于參數化建模思想提出針對彈性支承的分步優化設計方法[11]。栗勇等設計了一套高周疲勞試驗臺,對彈性支承高周疲勞性能進行試驗研究[12]。王永亮等使用試驗驗證鼠籠彈性支承疲勞故障機理并考慮加工方式對壽命的影響[13]。溫保崗等對籠條斷裂后鼠籠彈性支承的剛度特性進行研究[14]。王四季等采用主動控制方法,對鼠籠彈性支承下的轉子振動抑制進行了數值分析和試驗研究[15-17]。此外,雷沫枝等通過對彈性支承籠條動應力測試而獲得轉子的振動信號[18]。Inayat-Hussain等對帶有彈性支承的擠壓油膜阻尼器的轉子響應分叉問題進行研究[19]。Wang等針對彈性支承剛度提出了一種三維參數方程定義形變參數,并通過遺傳算法對其進行計算[20]。Ma等對帶有彈性支承橡膠阻尼器的轉子進行掃頻試驗,開展阻尼和剛度識別研究[21-22]。Degtiarev等對彈性支承剛度[23-24]和帶有彈性支承的擠壓油膜阻尼器進行參數識別試驗[25]。Kandhaswamy等針對航空發動機中彈性支承的剛度、強度和疲勞壽命的需求提出結構設計的方法[26]。

干摩擦阻尼器廣泛用于透平旋轉葉片減振,例如阻尼葉片和緣板阻尼器等。干摩擦包含界面黏滯-滑移等非線性特征,關于鼠籠彈性支承干摩擦阻尼器機理研究尚不充分。由于鼠籠彈性支承結構的材料參數和界面參數的不確定性,需要通過模態實驗進行參數識別和確定??紤]不同激振力和正壓力下的黏滯-滑移現象復雜,基于三維有限元模型的瞬態分析計算耗時嚴重。因此,有必要開展簡化梁單元模型的參數識別研究??紤]黏滯-滑移特性,本文提出一種雙線性彈簧梁單元簡化模型,模擬界面摩擦遲滯現象,兼顧計算精度和仿真效率。

本文建立了鼠籠彈性支承干摩擦阻尼試驗臺,開展干摩擦減振特性和參數識別研究。主要研究內容為:首先基于有限元模型和梁單元理論,開展彈性支承彎曲剛度分析,建立修正的剛度計算公式;建立鼠籠彈性支承干摩擦阻尼器減振試驗臺,改變接觸壓力,闡明最佳正壓力以及界面能量耗散特性;基于鼠籠彈性支承剛度特性及遲滯回線開展鼠籠接觸界面參數識別,提出了兩步參數識別方法,通過建立的雙線性彈簧模型對掃頻試驗激振力-位移響應遲滯回線進行驗證,與試驗結果吻合良好。

1 鼠籠剛度分析

鼠籠彈性支承是同心型擠壓油膜阻尼器的重要組成部分[27],近年來發生多起航空發動機阻尼環磨損和斷裂事故,造成巨大經濟損失[28]。鼠籠剛度設計不僅影響擠壓油膜阻尼器動力學性能和轉子響應,而且對結構疲勞壽命影響顯著。因此,有必要開展鼠籠彈性支承剛度分析。

1.1 有限元模型

首先建立鼠籠彈性支承有限元模型,使用Hypermesh軟件對鼠籠彈性支承進行網格劃分,鼠籠材料為40Cr,密度為7 850 kg·m-3,彈性模量為2.1×105MPa。

鼠籠是一種典型的彈性結構,如圖1所示,法蘭端固定的鼠籠可等效為懸臂桿梁結構。

L、b、h—籠條長度、寬度、厚度;F—橫向載荷。圖1 鼠籠結構示意圖Fig.1 Structure of the squirrel cage

因此,采用有限元方法進行靜力學分析,鼠籠彈性支承彎曲剛度k可表示為

(1)

式中:F為鼠籠自由端部施加的橫向載荷;δ為彈性支承響應位移。

改變鼠籠肋條寬度,彈性支承的結構參數見表1,圖2為鼠籠彈性支承的結構化網格。

表1 鼠籠彈性支承結構參數

圖2 鼠籠彈性支承的結構化網格Fig.2 Meshes of squirrel cage elastic support

分別計算不同鼠籠肋條寬度b=4, 6, 8, 10, 12 mm下的彎曲剛度,在法蘭側和端面分別添加全約束和載荷F,開展靜力學分析,根據式(1)得到的彎曲剛度隨寬度變化如圖3所示,可以看出剛度隨籠肋條寬度的增加而逐漸增大。以下剛度分析以有限元結果為參考值。

圖3 鼠籠彈性支承肋條寬度與彈性支承彎曲剛度的關系Fig.3 The relationship between the rib width of the squirrel cage and the bending stiffness of the elastic support

1.2 梁單元等效模型

將鼠籠等效為梁單元,文獻[29]采用彎曲剛度經驗公式為

(2)

式中:N為鼠籠條數;E為材料的楊氏模量。

此外,文獻[7]提出一種修正的理論公式

(3)

可以看出,式(2)和(3)均采用矩形截面梁單元進行簡化和修正,但是當鼠籠寬度較大時,采用矩形截面的誤差較大。將鼠籠等效為弧形截面梁,基于極坐標法得到了修正的鼠籠彎曲剛度計算公式。

以鼠籠中心為原點建立坐標系,鼠籠彈性支承截面如圖4所示。外力F作用方向為y軸負方向,假設第一個籠條邊與y軸正方向之間的夾角為θ,同一根籠條兩邊之間的夾角為α,兩個相鄰籠條之間的夾角為2π/N。

圖4 鼠籠彈性支承截面示意圖Fig.4 Cross section of squirrel cage elastic support section

每個籠條的截面為扇環形,y軸正方向開始順時針第i個籠條相對于x軸的慣性矩Ixi為

(4)

籠條兩邊之間的夾角α可以用如下公式估算

(5)

每根籠條的形心縱坐標為

(6)

使用移軸公式求出每根籠條相對于過其形心與x方向平行的軸的慣性矩如下式

(7)

整個鼠籠的剛度為所有籠條剛度之和

(8)

特別地,當N為4的倍數時

(9)

此時的鼠籠剛度各項同性,與外力作用方向無關。表2匯總了不同籠條寬度下的彈性支承剛度結果,其中λ為無量綱數鼠籠條寬度厚度比,λ=b/h;kFEM為仿真計算剛度大小;k1、k2和k3分別為式(2)、(3)和式(9)的計算結果。Δ1、Δ2和Δ3分別為3個理論公式和有限元分析的相對誤差。可以看出,彈性支承的剛度隨著鼠籠條寬度的增加而增加,但是理論公式的偏差也越來越大,這說明隨著籠條寬度的增加,不能簡單等效為梁單元的并聯。本文推導的理論公式的偏差Δ3相較于已有文獻公式偏差Δ1和Δ2更小。

進一步以λ作為變量進行分析,鼠籠彈性支承修正公式(9)與經典公式(2)(3)關于仿真計算剛度誤差對比如圖5所示。當λ=1.5時,本文所推導的剛度理論計算公式相較于仿真計算結果的誤差小于2%,而文獻[29]和[7]中公式的相對誤差Δ1、Δ2分別為5.24%和6.42%。

圖5 鼠籠彈性支承修正公式與經典公式誤差對比Fig.5 The error comparison chart of the corrected formula and the classic formula for the elastic support of the squirrel cage

2 干摩擦阻尼試驗流程

鼠籠彈性支承是航空發動機的重要非旋轉部件,干摩擦阻尼器可用于減少在啟停過程中轉子支承系統共振時的振動。本文在鼠籠端面添加干摩擦片并考慮不同正壓力時,進行鼠籠彈性支承的位移響應試驗和結果分析。

圖6為本試驗臺測試結構原理,主要包括輸入部分和試驗裝置本體部分,輸入部分包括信號發生器、功率放大器、激振器;試驗裝置本體包括重塊、靜摩擦片外框架、靜摩擦片、動摩擦片、鼠籠彈性支承、彈性支承外框架。此外還需要力傳感器、位移傳感器、數據采集卡測量和輸出信號。

圖6 鼠籠彈性支承試驗臺原理Fig.6 Schematic diagram of the squirrel cage elastic support test bench

鼠籠彈性支承干摩擦試驗流程如圖7所示,輸入部分主要通過信號發生器控制振動的頻率,通過信號發生器和功率放大器的共同作用控制振幅;在激振力作用下,彈性支承和動摩擦片一起振動,靜摩擦片安裝于靜摩擦片外框架內,由外框架控制靜摩擦片徑向位移,但靜摩擦片可以軸向滑動,通過重塊產生的重力提供摩擦副之間(靜摩擦片和動摩擦片)的壓緊力。通過改變激振力的大小和頻率,可以測量鼠籠在共振頻率附近的響應,從而得到干摩擦對于共振的影響規律;力傳感器和位移傳感器分別測量激振力和鼠籠響應,并通過NI9205數據采集卡進行采集和數據處理,得到不同接觸壓力下的頻響函數。

圖7 鼠籠彈性支承干摩擦試驗流程圖Fig.7 Flowchart of the squirrel cage elastic support dry friction test

圖8和圖9分別為試驗臺空載和添加重塊模擬干摩擦效應的實物圖。

圖8 鼠籠彈性支承實物圖(無靜摩擦盤)Fig.8 Squirrel cage elastic support without the rubbing disk

圖9 鼠籠彈性支承實物圖(有摩擦盤)Fig.9 Squirrel cage elastic support with the rubbing disk

3 鼠籠彈性支承干摩擦減振試驗

3.1 鼠籠彈性支承力-位移遲滯回線測量

干摩擦阻尼的阻尼由摩擦力所提供,由于摩擦力的存在,摩擦面的摩擦狀態有黏滯和滑移兩種狀態,這兩種不同的接觸面狀態也影響著鼠籠彈性支承的響應曲線從而影響減振特性,為研究彈性支承干摩擦減振特性,對鼠籠遲滯回線進行分析。

開展準靜態黏滯-滑移測量試驗研究,圖10是鼠籠彈性支承在5 Hz激振下的3組試驗結果。其中只改變正壓力負載的大小,當正壓力增大時,遲滯回線的橫向位移逐漸減小,即在輸入的激振力基本一致的情況下,彈性支承的位移響應逐漸減小,且在整個變化過程中,遲滯回線的四邊基本處于平行狀態,說明在鼠籠彈性支承激振狀態下,其運動的剛度基本保持不變。遲滯回線運動方向為1—2—3—4—1,其中點1為彈性支承響應位移負向最大點;點2為黏滯狀態轉變為滑動狀態轉變點;點3為響應位移正向最大點;點4為黏滯狀態轉變為滑動狀態轉變點。鼠籠彈性支承干摩擦阻尼器的運動狀態為:1—2為黏滯狀態;2—3為滑移狀態;3—4為黏滯狀態;4—1為滑移狀態。

圖10 5 Hz激振時鼠籠彈性支承改變正壓力試驗的力-位移遲滯回線Fig.10 The force-displacement hysteresis loop of a squirrel cage elastic support with 5 Hz excitation test

鼠籠彈性支承的響應力-位移曲線呈平行四邊形,存在遲滯現象,此時可以將此曲線簡化為一個雙線性剛度模型,如圖11所示,其中,當摩擦片處于黏滯狀態時,鼠籠彈性支承摩擦阻尼器的剛度為ks+km,當處于滑移狀態時剛度為km,兩條平行曲線封閉組合而成一個完整的平行四邊形。

圖11 鼠籠彈性支承簡化力-位移遲滯回線Fig.11 Simplified force-displacement hysteresis loop for squirrel cage elastic support

根據上述簡化模型對鼠籠彈性支承為5 Hz時改變正壓力試驗的結果進行參數識別,得到結果見表3。

表3 鼠籠彈性支承改變正壓力試驗雙線性剛度參數識別結果

3.2 最佳正壓力測量

彈性支承干摩擦阻尼器的阻尼主要由干摩擦提供,根據文獻[10]敘述,在改變彈性支承摩擦片正壓力時存在最佳正壓力,使得彈性支承減振效果最好,即系統等效黏性阻尼比最大,現對彈性支承進行掃頻試驗并改變彈性支承正壓力對最佳正壓力點附近進行試驗研究。

不同正壓力時鼠籠彈性支承頻響函數曲線如圖12所示。掃頻試驗在彈性支承一階共振點f為135 Hz附近進行,由圖12可以看出,當彈性支承在保證激振器信號發生器及功率放大器不變的情況下,改變彈性支承所受正壓力大小時,隨著正壓力的不斷增大,彈性支承的摩擦狀態由動摩擦狀態轉變為黏滯狀態,其中正壓力為14.21 N和20.58 N時為彈性支承的動摩擦狀態,24.50 N、27.93 N、30.87 N和34.30 N時為彈性支承動摩擦與黏滯轉變的中間狀態,最后兩組正壓力為37.24 N與40.67 N時彈性支承的摩擦狀態為黏滯狀態。

圖12 不同正壓力時鼠籠彈性支承頻響函數曲線Fig.12 Frequency response function curves of squirrel cage elastic support under different pressure

由于處于黏滯狀態無法使用半功率法計算阻尼比大小,現將非線性的摩擦阻尼等效為黏滯阻尼,即假設黏性阻尼ce,其在一個周期內做的功等于摩擦力所做的功,表達式為

(10)

式中:ΔE為一個周期所消耗能量;ω為激振器頻率;B為振動幅值。

(11)

式中:ζe為等效黏性阻尼比;ccr為臨界阻尼系數;m為彈性支承質量;ωn為彈性支承共振頻率。

聯立式(10)、(11)可得

(12)

計算得出鼠籠彈性支承改變正壓力試驗的等效黏性阻尼比見表4,其中Δe為等效黏性阻尼比相對空載增幅;σc為動摩擦片與靜摩擦片的接觸壓應力。

表4 鼠籠彈性支承改變正壓力試驗的等效黏性阻尼比對比

根據表4內容繪制出如圖13所示鼠籠彈性支承最佳正壓力曲線,由圖13中可以看出,隨著正壓力增大,彈性支承干摩擦阻尼器的等效黏性阻尼比會有一個先增大再減小的過程,即存在最佳正壓力,使得彈性支承干摩擦阻尼器的等效黏性阻尼比為最大值。

圖13 鼠籠彈性支承最佳正壓力曲線Fig.13 Optimum normal pressure curve for elastic support of squirrel cage

4 鼠籠彈性支承參數識別

4.1 鼠籠一維模型修正

由于三維模型中摩擦的瞬態計算耗時較長,通過3.1節中對鼠籠彈性支承的簡化和摩擦行為的分析,現由兩步修正將鼠籠彈性支承干摩擦阻尼器簡化為一維模型,圖14為模型修正的兩步流程。鼠籠彈性支承干摩擦阻尼器可以分為兩個部分,第一部分為鼠籠彈性支承本體;第二部分為摩擦盤部分。對于鼠籠彈性支承本體,可以將其簡化為一維梁單元;對于摩擦界面,可以簡化為雙線性彈簧模型,將非線性部分簡化為線性模型。

圖14 鼠籠彈性支承模型修正流程Fig.14 Model updating flow chart of squirrel cage elastic support model

鼠籠彈性支承干摩擦阻尼器的兩步簡化修正具體流程如下:首先對鼠籠彈性支承進行簡化,分析鼠籠彈性支承空載條件下的掃頻結果,計算頻響函數,提取10個峰值點,比較計算值與試驗值差值,采用優化算法對梁單元模型進行迭代計算,取計算中最佳優化點,根據最佳優化點建立梁單元模型;再對掃頻試驗中遲滯曲線進行參數識別,建立雙線性彈簧模型,將梁單元模型與雙線性彈簧模型結合即為鼠籠彈性支承干摩擦阻尼器的簡化模型。

第一步修正為彈性支承的一維梁單元模型簡化。優化算法利用響應面優化,設定梁單元長度H、梁截面半徑R、梁材料的楊氏模量E以及材料密度ρ4個參數,設定4個變量的范圍,計算400個計算點。再利用響應面優化中的多目標遺傳算法(MOGA),把試驗頻響函數的10個峰值點作為參考點,設置均方根誤差RMSE為0.02進行計算。圖15是計算的修正收斂曲線。

圖15 梁單元模型修正收斂曲線Fig.15 Beam element model correction convergence curve

(13)

選取其中誤差為0.01%的優化點,此時各參數值見表5。

對空載下優化修正參數進行驗證,使用修正簡化模型的掃頻頻響函數與實際試驗結果進行對比,圖16為空載時鼠籠彈性支承梁單元模型與試驗結果頻響函數對比,由圖16中可以看出,梁單元模型很好地擬合了彈性支承試驗的頻響函數曲線,均方根誤差為0.017。

圖16 空載時鼠籠彈性支承梁單元模型試驗與仿真頻響函數對比Fig.16 Comparison chart of frequency response function of beam element model test under elastic support of squirrel cage

第二步修正為鼠籠彈性支承干摩擦阻尼器摩擦界面的參數識別。3.1節中提到,鼠籠彈性支承阻尼器的遲滯回線存在雙線性的現象,可以對鼠籠彈性支承干摩擦阻尼器遲滯回線進行參數識別,提取雙線性剛度值,將摩擦界面簡化為雙線性彈簧單元,將鼠籠彈性支承簡化為梁單元,得到如圖17所示鼠籠彈性支承干摩擦阻尼器簡化模型。其中K0為梁單元剛度;K1為雙線性彈簧單元中的常值剛度;K2為雙線性彈簧單元中的摩擦剛度。激振時,當K1的彈性力小于Fslide時,簡化模型整體剛度為K0+K1+K2,即圖11中ks+km;當K1彈性力大于Fslide,簡化模型整體剛度變為K0+K2,即圖11中km。至此,兩步參數識別完成,鼠籠彈性支承干摩擦阻尼器三維模型可以簡化為一維帶彈簧的梁單元模型。

圖17 鼠籠彈性支承干摩擦阻尼器簡化模型Fig.17 Simplified model of dry friction damper with squirrel cage elastic support

4.2 計算結果對比

根據表5梁單元結構參數建立梁單元模型,開展了諧響應分析,其中阻尼比選取表4中14.21 N與24.50 N對應ξe,得到圖18所示滑移與黏滯狀態下頻響曲線對比圖。對其進行誤差分析,其中正壓力為14.21 N時共振頻率相對誤差Δf為0.72%,峰值相對誤差ΔB為0.78%。

圖18 滑移與黏滯狀態下試驗-仿真頻響函數對比Fig.18 Test-calculation frequency response function comparison chart

當對激振頻率為5 Hz時鼠籠彈性支承位移響應進行瞬態計算,時間步設置為0.001 s,網格單元尺寸設置為10 mm,三維模型計算0~1 s鼠籠彈性支承瞬態響應結果的時間t0為822 s,一維梁單元模型的計算時間t1為93 s,鼠籠彈性支承干摩擦阻尼器梁單元模型仿真計算時間比三維模型仿真計算時間相對減少了92.34%,計算效率顯著提高。

將瞬態仿真結果與試驗結果進行對比,曲線如圖19所示??梢钥闯?使用雙線性剛度簡支梁模型計算得出遲滯回線與試驗結果斜率基本一致,說明此模型能夠很好地模擬鼠籠彈性支承的摩擦行為。

(a)正壓力為14.21 N

(b)正壓力為20.58 N

仿真曲線與試驗曲線面積進行對比,正壓力為14.21 N時試驗與仿真誤差為8.76%,正壓力為20.58 N時試驗與仿真誤差為9.42%,兩種正壓力時試驗與仿真曲線的誤差均在10%以內。

5 結 論

本文開展了鼠籠彈性支承干摩擦減振試驗研究和參數識別,得到的結論如下。

(1)開展變結構參數的鼠籠彎曲剛度分析作為參考值,并基于極坐標系下的梁彎曲模型,建立了鼠籠彎曲剛度理論公式,相比已有理論公式精度更高。當寬厚比λ為1.5時,修正公式相較于經典公式(1)誤差減小了65.27%,相較于經典公式(2)誤差減小了71.65%。

(2)開展鼠籠彈性支承進行干摩擦減振試驗,得到了鼠籠干摩擦阻尼器的黏滯-滑移曲線。改變鼠籠彈性支承正壓力負載大小,通過等效黏性阻尼變化規律確定鼠籠彈性支承干摩擦阻尼器存在最佳正壓力,使得等效黏性阻尼比最大,即減振效果最好。

(3)建立雙線性彈簧模型,采用遺傳算法識別鼠籠彈性支承剛度;對鼠籠彈性支承干摩擦阻尼器黏滯-滑移曲線進行參數識別,得到鼠籠彈性支承干摩擦阻尼器摩擦界面特性。

(4)對參數識別結果進行試驗驗證和誤差分析,鼠籠彈性支承干摩擦減振器掃頻試驗誤差小于1%,遲滯曲線面積誤差小于10%。研究結果表明,雙線性梁單元模型能夠反映界面摩擦特性,且相比于三維有限元模型,計算時間減少了92.34%,計算效率大大提高,為鼠籠彈性支承干摩擦阻尼器的設計提供了一種理論計算模型。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲男人天堂2018| 91在线一9|永久视频在线| 国内嫩模私拍精品视频| 日韩欧美91| 亚洲天堂视频在线观看免费| 亚洲精品无码av中文字幕| 午夜老司机永久免费看片| 美女视频黄频a免费高清不卡| 国产99热| 日韩第一页在线| 亚洲视频在线观看免费视频| 国产 日韩 欧美 第二页| 欧美日韩一区二区在线播放| 国产成人夜色91| 最新日韩AV网址在线观看| 国产传媒一区二区三区四区五区| 精品国产一区91在线| AⅤ色综合久久天堂AV色综合 | 亚洲第一视频区| 国产激情无码一区二区APP | 久久这里只精品国产99热8| 亚洲日本一本dvd高清| 97超碰精品成人国产| 国产亚洲精品在天天在线麻豆| 无码粉嫩虎白一线天在线观看| 伊人久久综在合线亚洲91| 午夜福利在线观看入口| 国产激情国语对白普通话| 中文无码精品A∨在线观看不卡| 99久久精品免费看国产电影| 国产波多野结衣中文在线播放| 成人日韩精品| 国产主播喷水| 伊人精品成人久久综合| 精品久久高清| 国产男人天堂| 国内毛片视频| 在线观看国产精品第一区免费| 欧美天堂在线| 华人在线亚洲欧美精品| 国产91麻豆免费观看| 精品无码专区亚洲| 97在线公开视频| 欧美第二区| 91精品在线视频观看| 麻豆精品在线视频| 亚洲国产一区在线观看| 日韩在线永久免费播放| 国产三级毛片| 亚洲国语自产一区第二页| 久久精品一卡日本电影| 国产呦视频免费视频在线观看| 国产午夜一级淫片| 久久夜色精品| 色屁屁一区二区三区视频国产| 午夜不卡视频| 亚洲国产成人久久77| 国产91色在线| 国产一级在线播放| 久久精品日日躁夜夜躁欧美| 欧美日韩精品在线播放| 色妞永久免费视频| 72种姿势欧美久久久久大黄蕉| 日韩无码视频专区| 色悠久久久久久久综合网伊人| 91系列在线观看| 国产拍揄自揄精品视频网站| 欧美一区国产| 久久免费视频播放| 国产熟睡乱子伦视频网站| 毛片基地视频| 国产精品人莉莉成在线播放| 日韩精品久久无码中文字幕色欲| 国产女人在线视频| 99热国产这里只有精品无卡顿"| 亚洲综合网在线观看| 欧美一级夜夜爽www| 免费毛片a| 国产精品一区在线麻豆| 手机成人午夜在线视频| 国产第三区| 一区二区午夜|