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

航空發動機雙轉子試驗臺臨界轉速匹配優化

2021-10-20 10:30:58王東強陳坤旭李明明于賀春
機械設計與制造 2021年10期
關鍵詞:發動機優化模型

王東強,陳坤旭,李明明,于賀春

(中原工學院機電學院,河南 鄭州451191)

1 引言

試驗研究作為航空發動機動力學特性研究的重要手段之一,它對于航空發動機的設計、制造等技術研究至關重要,考慮到航空發動機試驗研究的高風險和高費用,目前模擬試驗臺便成為了試驗研究的常用方法[1]。隨著推重比的提高,航空發動機多為雙轉子結構,高低壓轉子間以及轉子系統與機匣的耦合影響日益增加,構成了較復雜的動力學特性[2]。因此對以航空發動機動力學特性為研究對象的試驗臺架而言,進行一定程度的簡化非常有必要,但簡化必然會帶來各種誤差,為了更真實的反映航空發動機的動力學特性,需保證基于簡化模型建立的試驗臺動力學特性與航空發動機動力學特性基本一致,其對航空發動機轉子動力學實驗領域的研究工作具有重要意義,也為后續的理論分析和實驗規律的探索建立了基準。對雙轉子系統進行臨界轉速的分析,是進行雙轉子系統動力學分析的基礎,臨界轉速的計算正確與否,直接影響后期雙轉子系統的瞬態響應、穩態響應、不平衡響應、穩定性分析、碰摩振動特性,以及不對中等動態特性的分析,因此,急需一種簡便、快速計算航空發動機雙轉子試驗臺與航空發動機原型機臨界轉速相匹配的優化算法,以便可以更好地對航空發動機的動力學特性進行試驗研究,同時也可縮短雙轉子試驗臺的設計周期。

從已公開的文獻來看,目前國內外,現有的航空發動機試驗臺大多是單轉子試驗臺,雙轉子試驗臺較少。此外,與國內的轉子試驗臺相比,國外的轉子試驗臺超臨界運行的工作轉速比國內的高很多。

同時,研究雙轉子試驗臺與航空發動機臨界轉速匹配性優化的文獻也較少。如文獻單轉子試驗臺主要對航空發動機低壓軸進行優化研究,轉子的最高轉速達到7900rpm[3];文獻單轉子試驗臺按Rolls-Royce RB401航空發動機進行設計,轉子的最高轉速為6000rpm[4-5];文獻建立的單轉子試驗臺主要通過調節轉子的彈性支撐從而研究轉子的臨界轉速[6-7];文獻采用了能夠滿足試驗要求的單輪盤、兩端支撐的單轉子試驗臺,轉子的臨界轉速主要是通過更換不同的輪盤來調節,額定轉速為12000rpm[8];文獻利用結構比較簡單的單輪盤、兩端支撐的單轉子試驗臺研究轉子的碰摩故障,速度較低[9];文獻研究中心的試驗臺為小型轉子試驗臺,主要是用于氣箔軸承動力學特性研究,運行速度為20000rpm[10-11];文獻[12]和文獻[13]建立的試驗臺為雙轉子試驗臺,其能達到的最高轉速都較小;文獻[14]的雙轉子試驗臺較小,長度不超過1m;柏林工業大學的雙轉子試驗臺的高低壓轉子的工作轉速相對較高,分別為12500rpm和6000rpm[15-16]。

文獻應用ANSYS有限元軟件計算了雙轉子系統的臨界轉速以及振型[17];文獻建立了雙轉子軸承系統的動力學模型,研究了不同轉速比情況下雙轉子系統的臨界轉速[18];文獻以轉子-軸承系統的不平衡響應幅值為優化目標,優化了該系統的臨界位移,研究對象為單轉子系統[19];文獻以多階臨界轉速及動力響應為約束變量對一單轉子系統的多階振型進行了優化[20];文獻運用混合遺傳算法對一轉子系統的臨界轉速進行了優化設計,優化對象為單轉子系統[21]。

文獻以大型汽輪發電機組為研究對象,計算了其臨界轉速,研究了單轉子系統的優化設計[22];文獻以球軸承-轉子系統為研究對象,對該系統的固有頻率與臨界轉速的分析及其結構的優化進行了研究[23];文獻在多支撐轉子-軸承系統試驗臺設計的過程中,為了更好的反映多支撐軸承-轉子系統的動力學特性,盡量保證試驗臺的前三跨轉子前兩階臨界轉速與實際轉子一致,研究對象仍為單轉子系統[24];文獻在風扇轉子試驗器的設計中,逐個分析了轉軸內外徑、輪盤外徑與質量、軸承的剛度對試驗器臨界轉速的影響,通過調整這些參數實現試驗器臨界轉速等特性的控制[25];文獻提出了一種雙轉子實驗臺結構,逐個分析了支承剛度、轉速比、輪盤的極轉動慣量、長徑比等因素對臺架臨界轉速的影響,并據此對實驗臺架的臨界轉速做了調整[26];文獻逐個分析了轉子的支撐方式、支點跨距、轉軸直徑大小、轉子質量分布以及支撐剛度對試驗器轉子的臨界轉速的影響,據此對試驗臺轉子各階臨界轉速進行了調節[1]。

綜合上述文獻,國內外高校和研究院所搭建的航空發動機實驗平臺大部分都只遵循結構相似原則,對于模擬汽輪機和渦輪的輪盤在轉子上的定位比較隨意,計算得到的雙轉子臨界轉速與航空發動機臨界轉速的實測值存在著較大的誤差,雙轉子試驗臺的尺寸普遍較小,轉速相對較低;雙轉子系統臨界轉速優化過程中多單獨根據一組設計參數進行逐步優化,按照該方法再分析某個參數對其臨界轉速的影響,進行調節時會對多階臨界轉速都會施加影響,逐個參數的反復調節必然會消耗大量的時間;另外現代航空發動機向高速化和柔性化方向發展,尤其是中小型航空發動機,更多地使用薄壁板殼結構,具有復雜而密集的振動模態,轉子工作轉速往往位于其幾階臨界轉速以上,單純通過經驗和試湊進行多階臨界轉速匹配性優化,工作量過大,可靠性也不易保證。

鑒于此,提出了雙轉子試驗臺和航空發動機原型臨界轉速匹配性優化方法,研究對象為某型號航空發動機雙轉子,低壓轉子額定工作轉速為12000rpm,高壓轉子額定工作轉速為15000rpm,根據其支撐方式和結構形式,搭建了其轉子-軸承-機匣耦合系統簡化模型,基于有限元方法求解了模型臨界轉速,運用變換哈墨斯利算法實現模型與原型機臨界轉速匹配性優化。

2 雙轉子航空發動機整體動力學模型

2.1 轉子模型

根據某型號航空發動機的支撐方式和結構形式等建立轉子-軸承-機匣耦合系統簡化模型,如圖1所示。該系統的支撐結構為五點支撐,軸承支撐的編號從左到右依次是2,6,10,13,編號10為中介軸承,軸承選用航空發動機同類型軸承,編號2處軸承類型為深溝球軸承,編號6和10處軸承類型為角接觸球軸承,編號10和13處軸承為圓柱滾子軸承,系統由低壓轉子、低壓壓氣機模擬盤、高壓轉子(空心)、高壓壓氣機模擬盤、高壓渦輪模擬盤及動力渦輪、低壓渦輪模擬盤和機匣組成、高壓轉子和低壓轉子靠中介軸承連接。

圖1 轉子-軸承-機匣耦合系統簡化模型Fig.1 Simplified Model of Rotor-Bearing-Casing Coupling System

為了使得模型的復雜程度減小,方程的求解更高效,模型的建立過程中遵循以下假設:

(1)低壓轉子和高壓轉子由于都具有較大的長徑比,因此可用Timoshenko梁單元表示,只考慮這些轉子軸的橫向彎曲振動,不考慮轉子軸的扭轉振動和軸向振動。

(2)假設轉盤以及轉子軸的材料是均勻連續和各向同性的,并在彈性范圍內服從Hook定律。

(3)低壓壓氣機、高壓壓氣機、高壓渦輪、動力渦輪以及低壓渦輪都以模擬盤的形式表示,以集中質量和轉動慣量的形式作用于模擬盤的質心上,考慮盤的陀螺效應。

(4)假設在工作速度范圍內機匣和轉子軸產生彎曲耦合,其橫截面的形狀為圓形,故可將機匣按照固定的梁單元表示。

(5)軸承被簡化為彈簧-阻尼系統,不考慮軸承的交叉剛度。

2.2 模型運動方程

從圖1可知,模型的絕對坐標系為o-xyz,坐標軸o-z與高壓轉子和低壓轉子的軸線重合,建立的轉子-軸承-機匣耦合系統簡化模型主要包括模擬盤、柔性軸、機匣和支撐四類部件,根據文獻[27]、[28]按照部件類型完成單元動力學方程的建立,根據力學關系進行部件組裝,最終得到轉子-軸承-機匣系統整體運動微分方程。

2.2.1 模擬盤

模擬盤假設為剛性圓盤,不考慮其應變能,僅考慮其動能,取模擬盤的橫向線位移為x,y和夾角θx,θy為廣義坐標。設盤的質量為md,直徑轉動慣量為Jd,極轉動慣量為Jp,求取模擬盤的動能,運用Lagrange方程得到模擬盤的運動方程為:

式中:Md-模擬盤的質量矩陣,J-陀螺矩陣,模擬盤軸心的廣義位移向量為q1d=[x,θy]T,q2d=[y,-θx]T,Q1d和Q2d為相應節點處的廣義力向量包括不平衡力、圓盤相鄰軸段的力和力矩以及支撐的約束力等。則上述矩陣可以表示為:

2.2.2 軸段微分方程

建立軸段單元如圖2所示,取軸段AB兩端的節點的位移為廣義坐標即:

圖2 彈性軸段Fig.2 Elastic Shaft Segment

通過Lagrange方程可以得到軸段AB在給定坐標系下的運動微分方程為:

式中:MST-軸段的單元移動慣性矩陣;MSR-單元轉動慣性矩陣;Ks-單元剛度矩陣;Js-單元回轉矩陣;Q1s和Q2s-廣義力向量;設軸段的長度l,軸截面的面積為A,軸的半徑為r,密度為ρ,軸段的單位長度的質量為α即α=ρA,則上述單元矩陣可以表示為:

2.2.3 軸承運動微分方程

如圖3所示,假設軸承座中心的節點處的坐標為xs,ys,對應的軸徑中心節點的坐標為xj,yj,則軸承座的運動方程為:

圖3 軸承模型Fig.3 Bearing Model

2.2.4 機匣運動微分方程

機匣運用固定的梁單元進行建模,考慮了其剪切效應和轉動慣量,采用有限元法可以得出機匣的運動方程為:

2.2.5 系統運動微分方程

按照高壓轉子、低壓轉子、模擬盤、軸承以及機匣之間耦合節點的對應關系,以及界面上的力的相互作用關系,從而合并各類部件的運動方程,便可以得出雙轉子系統的運動微分方程為:

式中:M1、M2-高低壓轉子系統的總質量矩陣;K1、K2-高低壓轉子系統的總剛度矩陣;Ω1J1、Ω2J2-高低壓轉子系統的總陀螺矩陣;Q1,Q2,Q3,Q4-外部激勵;h、l-高低壓轉子的節點數,Ω1和Ω2-高低壓轉子的轉速,令λ=Ω2/Ω1即低高壓轉子轉速比,K12-中介軸承10所產生的耦合剛度可表達為:

式中:-kxx、-kyy-(4q-3,4p-3)和(4q-2,4p-2)處;p和q-中介軸承節點對應的高低壓轉子節點數。

則,系統運動微分方程可以表示為:

令,Qu=0,Qv=0則可以得到系統微分方程對應的齊次方程的解為:

將方程(15)代入方程(14)可以得出:

從而得到頻率方程為:

根據式(17)便可以求出當Ω1=ω時,轉子-軸承-機匣耦合系統的臨界轉速,此時得出了高壓轉子為主激勵時系統臨界轉速。

3 轉子系統臨界轉速計算

表2 高壓轉子軸參數表Tab.2 Parameters of High Pressure Rotor Shaft

表3 模擬盤參數表Tab.3 Parameter Table of Analog Disk

表4 彈性支撐參數Tab.4 Elastic Support Parameters

表5 機匣計算參數Tab.5 Casing Calculation Parameters

結合圖1所示的轉子-軸承-機匣耦合系統模型,根據表1~5中的計算參數,編寫系統臨界轉速計算程序,高壓轉子與低壓轉子的轉速比為1.25,分別得到以低壓轉子為主激勵和高壓轉子為主激勵的雙轉子系統坎貝爾圖,如圖4、圖5所示。

表1 低壓轉子軸參數表Tab.1 Parameters of Low Pressure Rotor Shaft

圖4 低壓轉子為主激勵的坎貝爾圖Fig.4 Campbell Diagram with Low Pressure Rotor as Main Excitation

圖5 高壓轉子為主激勵的坎貝爾圖Fig.5 Campbell Diagram with High Pressure Rotor as Main Excitation

如圖4所示,表示了高低壓轉子以1.25的轉速比運轉,轉子整體系統的正反進動曲線隨低壓轉速變化圖,圖5則表示轉子整體系統的正反進動曲線隨高壓轉速變化圖。根據臨界轉速的定義,等速線與正進動曲線的交點即為臨界轉速值,得到模型由高低壓轉速激起的前四階臨界轉速計算結果,如表6所示。

表6 簡化模型臨界轉速計算結果Tab.6 The Critical Speed Calculation Results of the Simplified Model

4 航空發動機臨界轉速測試

模型的簡化必然會帶來誤差,為了對比模型臨界轉速與實測數據的誤差以及為后續臨界轉速優化建立基準,通過三維轉速譜圖分析法,對航空發動機整機的臨界轉速進行了測試。

4.1 測點分布與測試系統

測點分布如圖6所示,共放置5個傳感器,1和3測點分別為風扇前支點和中介機匣位置,放置磁電式速度傳感器,2、4和5測點分別為附件機匣位置、減速器和低壓渦輪支點,放置三向加速度傳感器,信號采集設備為美國DP信號采集系統,整機測試系統如圖7所示。

圖6 航空發動機整機測點布置Fig.6 Layout of Measuring Points for the Whole Aircraft Engine

圖7 航空發動機整機測試系統圖Fig.7 Test System Diagram of Aircraft Engine

4.2 測試結果及其與計算結果誤差分析

測試得到三維轉速譜圖,其橫坐標為頻率,縱坐標為以一定轉速相隔的響應功率譜,從圖中的功率譜曲線峰值走向以及變化情況便可以初步確定臨界轉速范圍,若峰值組成的“山脈”的振動頻率與轉速成正比關系,且在某一轉速范圍內其峰值明顯,則該轉速即為臨界轉速,根據該判斷方法,便得到整機臨界轉速,測試結果及其與模型計算結果的誤差,如表7所示。

表7 實測臨界轉速及其與模型計算結果誤差Tab.7 Measured Critical Speed and Its Error with Model Calculation Results

從表7可知簡化模型的臨界轉速與航空發動機實測數據高階臨界轉速相差較大,超過了允許誤差5%,為了滿足兩者動力學特性的相似性,需要對模型的臨界轉速進行優化。

5 轉子系統臨界轉速優化

在臨界轉速不滿足設計要求時,可改變轉子系統的結構參數,例如轉軸的長度、直徑、輪盤的尺寸等,或者軸承以及機匣的支承剛度,從而實現模型臨界轉速的調整[21]。為了更真實的反映航空發動機動力學特性,軸承選用航空發動機同類型軸承,不改變系統支撐剛度以及軸的直徑。根據表7中的數據,模型的高階臨界轉速需進行較大程度的調整,由于輪盤的轉動慣量所產生的陀螺效應對轉子的高階臨界轉速影響較大,因此優化輪盤的大小對于調節轉子的臨界轉速具有重要作用[1]。此外由于試驗臺場地的限制,轉子軸長度需滿足一定約束條件。考慮到采用人工調整的方式進行調節,無法保證優化的效率及準確性,采用變換哈墨斯利(Shifted Hammersley)算法對模型的臨界轉速進行優化。

5.1 變換哈默斯利算法構建原理

變換哈墨斯利算法優化方法是一種適用于所有類型樣本生成的一種抽樣方法,根據根逆函數構建而成,過程如下:

任何整數n都可以用下面的公式表示為一個數字序列n0,n1,n2,……,nm

如果一個整數用基數R表示,則式(18)可表示為:

逆根函數是指通過對式(19)中的關于小數點的數字順序進行倒轉,從而得到(0,1)中的分數的函數如下:

因此,對于k維搜索空間,Hammersley點由以下公式給出:

其中,i=0,…,N表示樣本點,為了使得Hammersley采樣器生成的點更靠近k維超立方體的原點,提出一個點移動過程,將所有Hammersley點按式(22)進行移動

5.2 臨界轉速優化設計模型

假設系統目標函數為分別以高低壓轉子為主激勵的臨界轉速,0<nh1<nl1<nh2<nl2<nh3<nl3<nh4<nl4,nhi和nli分別是以高低壓轉子為主激勵的第i階臨界轉速,約束變量為低壓轉子軸的長度,lil<li<liu(i=1,2,…,s),lil和liu分別為約束變量的上下限,ckl<ck<cku(k=1,2,…,t),ckl和cku分別為設計變量的上下限,則該優化問題可以表達為:式中:n*hi、n*

li-優化目標函數值即分別以高低壓轉子為主激勵的實測臨界轉速值。

5.3 靈敏度分析

優化的目標函數為系統中分別以高低壓轉子為主激勵的前四階臨界轉速,目標值為航空發動機真機的實測數據,設計變量為低壓壓氣機、高壓壓氣機、高壓渦輪、動力渦輪、低壓渦輪的轉動慣量,考慮到試驗臺場地的限制,將低壓轉子的軸段長度也同時作為設計變量,由于低壓轉子的軸段數較多,考慮到實際工況,排除左右兩端的聯軸器連接軸段,選擇如圖8所示的各軸段長度進行靈敏度分析。首先建立靈敏度分析的數學模型,設α為系統臨界轉速值,β1、β2、…、βm為轉子系統的各個參數值,則系統的特征方程為:

圖8 低壓轉子部分軸段長度Fig.8 Length of Shaft Segment of Low Pressure Rotor Part

對公式(25)兩邊微分可以得到:

由式(26)可得,系統的參數的增量是δβk(k=1,2,…,m),其臨界轉速的微小增量是δα,其轉子系統的特征行列式的改變量是δΔ。

令δβj≠0,δβk=0(k≠j,k=1,2,…,m),αl為轉子系統某一階的臨界轉速,則式(26)可表示為:

若δΔl=0,則轉子系統的靈敏度值為:

根據sl j的數值,便可以得出對轉子系統臨界轉速影響較大的參數。對圖8所示的各個軸段長度進行靈敏度分析,其結果如表8所示。

表8 低壓轉子軸段靈敏度分析結果Tab.8 Sensitivity Analysis Results of Low Pressure Rotor Shaft Segment

表中的系統靈敏度為無量綱,負號表示隨著系統參數增大,系統臨界轉速減小,軸段單位為m,臨界轉速單位為rpm。由于優化的目的是縮短低壓轉子軸段長度,減小高階臨界轉速,因此靈敏度選擇標準為較大的正值,綜合分析表中數據知N、K、S三段的靈敏度分析結果較符合選擇標準。

5.4 臨界轉速優化過程

系統優化流程,如圖9所示。根據初始數據建立了轉子-軸承-機匣耦合系統的有限元模型,求解系統的臨界轉速,優化進行之前設定了優化目標優先順序以及允許優化誤差,設定了約束變量以及設計變量范圍,優化過程中,根據優化算法進行設計變量的修改,再依次進行優化結果評估和判斷是否滿足約束條件的要求,若有一者不滿足,則重新修改設計變量,若兩者都滿足,則輸出解,若所有樣本數據都無法滿足,則重新定義設計變量的取值空間,重新進行優化計算,直到找到符合要求的解集,最后根據優化目標的優先順序得到最優解。

圖9 系統優化流程圖Fig.9 System Optimization Flow Chart

6 優化結果分析

優化完成后得出轉子分別以高低壓轉子為主激勵的坎貝爾圖如圖10、圖11所示。從而得到分別以高低壓轉子為主激勵的臨界轉速及其與實測數據的誤差,如表9所示。系統的設計參數的優化結果,如表10所示。

圖10 優化后低壓轉子為主激勵的坎貝爾圖Fig.10 Campbell Diagram of Optimized Low Pressure Rotor with Main Excitation

圖11 優化后高壓轉子為主激勵的坎貝爾圖Fig.11 Campbell Diagram of Optimized High Pressure Rotor with Main Excitation

表9 系統臨界轉速優化結果Tab.9 System Critical Speed Optimization Results

由表9和表10可得出系統臨界轉速優化結果與航空發動機真機實測數據相差較小,遠小于允許誤差,優化后低壓轉子的長度比優化前低壓轉子長度減少了9.7%,滿足了設計要求。從表中數據可以得出根據優化后參數進行試驗器的設計可以很好的反映航空發動機真機的動力學特性,從而為后續的理論分析和實驗規律的探索建立基準。

表10 系統設計參數優化結果Tab.10 Optimization Results of System Design Parameters

7 結論

基于航空發動機的結構,建立了轉子-軸承-機匣耦合系統的簡化模型,運用有限元方法求解了高低壓轉子轉速比為1.25時系統的臨界轉速,為了滿足計算結果與實測數據的匹配,結合變換哈默斯利算法進行了臨界轉速的優化。結果表明,雙轉子試驗臺的臨界轉速和航空發動機實測數據誤差在5%以內,優化后低壓轉子的長度比優化前低壓轉子長度減少了9.7%,滿足了設計要求,簡化模型可以較好的反映航空發動機真機的動力學特性,證明了該優化方法的有效性,對航空發動機試驗臺的設計具有指導意義。

猜你喜歡
發動機優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
發動機空中起動包線擴展試飛組織與實施
3D打印中的模型分割與打包
新一代MTU2000發動機系列
主站蜘蛛池模板: av午夜福利一片免费看| 2021精品国产自在现线看| 日本日韩欧美| 久久精品国产电影| 亚洲免费人成影院| 欧美日韩国产精品综合| 日韩欧美中文字幕在线精品| 思思热精品在线8| 亚洲成A人V欧美综合| 亚洲色成人www在线观看| 免费三A级毛片视频| 精品午夜国产福利观看| 欧美亚洲国产精品久久蜜芽| 一级福利视频| 国产成人AV男人的天堂| 日本精品αv中文字幕| 亚洲香蕉伊综合在人在线| 久久这里只有精品23| 都市激情亚洲综合久久| 免费在线a视频| 亚洲成人一区二区| 国产成人一级| 亚洲天堂网在线视频| 色香蕉影院| 欧美日韩中文国产| 日韩AV无码一区| 中文国产成人精品久久| 亚洲无码高清视频在线观看| 精品久久久久久久久久久| 日韩成人高清无码| 国产丰满大乳无码免费播放| 国产黄网永久免费| 亚洲欧美天堂网| 亚洲国产中文精品va在线播放 | 国产午夜无码片在线观看网站 | 小说 亚洲 无码 精品| 老司国产精品视频91| 日本在线国产| 97成人在线观看| 成人午夜视频免费看欧美| 国产精品99久久久久久董美香| 成人在线第一页| 亚洲第一香蕉视频| 最新国产精品第1页| 毛片网站观看| 欧美日韩激情在线| 不卡无码网| 91外围女在线观看| 欧美爱爱网| 国产欧美中文字幕| 亚洲永久视频| 欧美伦理一区| 欧美日本在线观看| 奇米影视狠狠精品7777| 午夜视频在线观看区二区| 国产十八禁在线观看免费| 久久久久中文字幕精品视频| 国产靠逼视频| 国产成本人片免费a∨短片| 精品人妻一区无码视频| 极品私人尤物在线精品首页 | 国产成年女人特黄特色毛片免| av在线人妻熟妇| 97人人做人人爽香蕉精品| 亚洲国产中文精品va在线播放| 久久久久久尹人网香蕉 | 欧美激情一区二区三区成人| 91啦中文字幕| 中文一区二区视频| 欧美成人国产| 伊人久久大线影院首页| 欧美成人午夜在线全部免费| 免费人成网站在线观看欧美| 黄色污网站在线观看| 又粗又大又爽又紧免费视频| 婷婷综合色| 老汉色老汉首页a亚洲| 香蕉视频国产精品人| 中文字幕不卡免费高清视频| 亚洲天堂网在线播放| 制服丝袜在线视频香蕉| 试看120秒男女啪啪免费|