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

基于整機瞬態熱應力場重型燃氣輪機轉子的壽命評估

2015-10-28 11:23:33歐文豪石清鑫
中國機械工程 2015年7期

歐文豪 袁 奇 石清鑫

西安交通大學,西安,710049

基于整機瞬態熱應力場重型燃氣輪機轉子的壽命評估

歐文豪袁奇石清鑫

西安交通大學,西安,710049

針對某F級重型燃氣輪機,根據機組的實際運行監測數據,提取得到該機組在冷態啟動、熱態啟動兩種工況下各測點溫度隨時間的變化曲線,分析得到轉子各處的換熱邊界條件,采用有限元方法計算該燃氣輪機轉子在不同啟動工況下的整機瞬態溫度場及應力場,并基于該瞬態熱應力分析結果對轉子的低周疲勞壽命損耗進行評估。計算結果表明:最大熱彈性應力值一般都出現在透平第一級輪盤的出氣側底部圓角處及壓氣機后三級輪盤的中間部位的圓角區域;氣流參數變化的快慢對轉子應力分布影響極大;冷態啟動一次的壽命損耗最高達0.02%,大于熱態啟動的對應值。研究結果為輪盤結構設計與機組運行規程的優化提供了重要參考。

燃氣輪機轉子;啟動工況;有限元方法;應力場;壽命評估

0 引言

燃氣輪機轉子多采用盤式拉桿轉子結構,工作在高溫、高轉速的惡劣環境下,高溫環境導致轉子內部存在不均勻的溫度分布,會存在大的熱應力,并在啟停過程中不斷變化,而轉子在額定工況穩態運行及各啟停過渡工況下的應力狀態與轉子的安全可靠性密切相關。針對轉子復雜的換熱邊界條件,國內外學者從理論、數值模擬、實驗方面對旋轉結構內的流動換熱現象進行了大量的研究。文獻[1]分析了蒸汽輪機疲勞與機組運行方式的關系,深入討論了如何通過熱應力監測來進行有效壽命管理。文獻[2-3]分析了燃機輪盤腔體內冷卻空氣的流動及換熱情況和輪盤溫度分布及熱應力分布的關系,重點討論了預旋系統和旋轉的外周流入和流出的冷卻空氣如何影響旋轉腔體附近區域乃至整個輪盤的溫度場。文獻[2-6]將汽輪機、燃氣輪機內的復雜流動結構進行一定的簡化,抽取典型的流動換熱結構特征來進行研究。在此基礎上,國內外學者在采用有限元數值模擬轉子熱應力及影響轉子熱應力因素等方面進行了大量的研究。文獻[7]以某電廠600 MW超臨界機組汽輪機轉子為研究對象,運用ANSYS有限元軟件建立二維模型,計算出不同平均溫升率下機組啟動時轉子熱應力值,與理論公式通過最小二乘法進行曲線擬合,從而確定無中心孔轉子的時間修正因子和形狀因子,然后根據低周疲勞曲線計算出轉子的壽命損耗。文獻[8]通過對轉子剖面溫度變化的實時監測,確定了汽輪機轉子的最大應力監測面,并對此監測面進行了分析。文獻[9]采用電廠實際啟停機曲線對600 MW汽輪機組的高中壓轉子分別進行了三維瞬態溫度場及非線性應力場分析,并利用30Cr1MoV材料的低周疲勞曲線對機組在實際啟停過程中的壽命損耗進行了估算。雖然關于轉子壽命評估和管理方法的研究已有很多,但其研究對象主要局限于蒸汽輪機,相較蒸汽輪機最高600 ℃的進氣溫度,F級燃氣輪機透平進口溫度高達1400 ℃左右,因此由溫度分布不均產生的熱應力也遠比蒸汽輪機嚴重。同時,為了降低高溫對轉子的不利影響,燃氣輪機轉子中均有冷卻通道與冷卻空氣,使得燃機透平結構與換熱邊界條件更為復雜。但由于其技術的保密性,對燃氣輪機轉子的熱應力及疲勞壽命進行研究的文獻很少。因此,采用有限元方法計算分析燃氣輪機轉子的熱應力并評估轉子的低周疲勞壽命損耗具有理論意義與實際工程價值。

本文以某F級重型燃氣輪機拉桿轉子為研究對象,建立了軸對稱轉子有限元模型;根據機組設計參數及相關監測數據計算得出轉子各個部位的換熱系數,得到符合實際的初始條件和邊界條件,采用ANSYS有限元軟件計算分析了某F級重型燃氣輪機轉子在冷態啟動和熱態啟動工況下的瞬態溫度場及應力場,獲得了兩種啟動工況下轉子的溫度場及應力場的變化規律;最后,基于瞬態分析結果對轉子的低周疲勞壽命損耗進行了評估。

1 研究模型及方法

1.1拉桿轉子有限元模型

某F級重型燃氣輪機主要由壓氣機、燃燒室和透平三個部分組成。壓氣機為17級、壓比為17的軸流壓氣機,用12根周向均勻分布的長拉桿連接;4級透平輪盤也用12根稍短的拉桿連接。

圖1 某F級重型燃氣輪機的結構示意圖

實際的轉子模型結構復雜,細微結構眾多,建模和邊界條件的簡化處理如下:忽略氣流力、重力、拉桿預緊力及扭矩等對轉子的強度影響,主要考慮轉子的熱應力及離心應力;對葉片及輪緣進行等效處理,將每一級轉子的葉片及圍帶等去掉,在葉片根部添加假想連續環狀體,與轉子作為一個整體;簡化對計算結果影響很小的圓角等細小結構,重要部位的幾何模型中的圓角和倒角均按照轉子精加工圖處理,以正確反映轉子的應力集中;整根轉子的左右端面和轉子中的封閉空腔結構做絕熱處理,外表面均作為已知換熱系數及周圍流體定性溫度的第三類邊界條件。

建立有限元計算模型,如圖2所示。由于轉子是盤式拉桿轉子,它是通過拉桿將輪盤組裝成一整體的,其實際轉子接觸面眾多,實際轉子拉桿預緊力足夠大,轉子的各級輪盤之間的接觸關系可以當作整體。

圖2 某F級重型燃氣輪機轉子的軸對稱有限元模型

1.2熱邊界條件的確定

由于轉子是高速旋轉的,依照現有的技術,還沒有有效的手段來實時監測運行中的燃氣輪機轉子與氣流之間的對流換熱狀況,一般都是通過相似性實驗對相似情況下的對流換熱系數進行經驗關系擬合。本文根據轉子結構的特點、氣流的流動特性及《航空發動機設計手冊》[10]的相關內容,對燃氣輪機轉子熱邊界條件進行了相應的簡化處理:①將葉片及輪緣進行等效處理,在輪盤外緣處添加與葉根寬度一致的不同密度材料的假想連續環狀體,在建模時將其加到輪盤的外緣處,與輪盤做成一體;②整根轉子左右端面做絕熱處理;③僅轉子內部空腔結構與冷卻空氣直接接觸的轉子表面給定換熱條件,其他內表面均給定絕熱邊界;④轉子的外表面均作為已知換熱系數及周圍流體定性溫度的第三類邊界條件。

根據簡化處理將燃機轉子的換熱邊界條件分為以下4個基本類型:輪盤側面換熱[10]、光軸處換熱[11]、壓氣機輪緣處換熱[11-12],透平輪緣處換熱。前三種換熱類型的準則式選取見表1。表1中,h為換熱系數,λc為氣流的導熱系數,Nu為努塞爾數,Re為雷諾數,u為外圓Rb處的圓周速度,ν為氣流運動黏度系數,Ra為光軸外半徑,R0為輪緣半徑,λ為葉片材料的導熱系數。由于透平的葉片結構復雜,且內部還有冷卻空氣,故其換熱系數很難確定。本文依據熱平衡原理,綜合考慮燃氣與葉片及葉片平臺的換熱、冷卻空氣對葉片的換熱、冷卻空氣對葉根間隙的換熱及葉片與輪緣的接觸熱阻等因素,采用迭代方法計算得到透平各輪緣處的等效換熱系數,其定性溫度取為燃氣的動葉進出口溫度的平均值。

表1 各換熱類型準則式

為了獲得燃氣輪機轉子的整體溫度場及應力場隨機組啟停的變化規律,本文進行了以下簡化處理:①轉子各表面換熱系數按隨轉速及流量變化而變化考慮;②基于比例系數,根據各啟停工況已有的監測點數據,推算得到各級特征面處的氣流溫度隨時間的變化規律。

1.3主要參考點

為便于直觀地分析結果,明確轉子危險點的應力隨啟停過程的變化規律,分別對溫度及應力較高的幾級輪盤定義了不同位置的參考點,用A,B,…,G表示,具體位置如圖3所示。

(1)點A位于透平第一級輪盤的輪緣處,其溫度代表透平第一級輪盤輪緣處溫度水平;

(2)點B位于透平第一級輪盤出氣側底部圓角處,該處應力集中明顯,是轉子壽命監視的危險點;

(3)點C位于透平第一級輪盤軸心線中點處,其溫度代表了透平第一級輪盤輪心處的溫度水平;

(4)點D位于壓氣機第十六級輪盤的外緣處,其溫度代表了壓氣機第十六級輪盤輪緣處的溫度水平;

(5)點E位于壓氣機第十六級輪盤的出氣側底部圓角處,該處有較大的應力集中,其應力變化代表壓氣機轉子應力隨啟停過程的變化規律;

(6)點F位于壓氣機第十六級輪盤的軸心線中點處,其溫度代表壓氣機第十六級輪盤輪心處的溫度水平;

(7)點G位于壓氣機第十五級輪盤的進氣側內腔圓角處,該處存在較大的應力集中,是壽命監視的危險點。

圖3 某F級重型燃氣輪機轉子各參考點具體位置示意圖

2 冷態啟動工況的計算結果與分析

冷態啟動時間大約210 min,但為了較好地模擬熱應力隨啟動過程的變化規律,取冷態啟動工況計算時間為334 min。根據已有的一組機組運行監測數據,提取了機組的冷態啟動工況下主要參數隨時間變化關系曲線,如圖4、圖5所示。

圖4 冷態啟動工況燃機轉速及功率隨時間變化曲線

1.壓氣機入口 2.壓氣機出口 3.透平出口 4.透平第二級靜葉腔室 5.透平第三級靜葉腔室 6.透平第四級靜葉腔室 7.透平第四級下端腔室 8.轉子冷卻空氣入口圖5 冷態啟動工況燃機主要監測點溫度隨時間變化曲線

2.1冷態啟動工況的計算結果

根據啟動過程的參數動態變化,將計算開始時刻的監測數據作為已知熱邊界條件,以停機72 h后的轉子溫度場分布作為初始條件,再通過ANSYS瞬態熱分析計算得到任意時刻轉子溫度場分布,綜合考慮離心力載荷與溫度載荷,計算得到任意時刻轉子綜合加載等效應力分布。如圖6~圖13所示。

圖6 初始時刻轉子溫度場分布

圖7 啟動后60 min轉子溫度場分布

圖8 啟動后334 min(計算終止時刻)轉子溫度場分布

1.外緣處D點 2.內徑處F點 3.兩點的溫差(a)壓氣機第十六級輪盤

1.外緣處A點 2.內徑處C點 3.兩點的溫差(b)透平第一級輪盤圖9 冷態啟動工況輪盤參考點溫度隨時間變化曲線

圖10 啟動后60 min轉子綜合加載等效應力分布

圖11 啟動后180 min轉子的綜合加載等效應力分布

圖12 啟動后334 min轉子的綜合加載等效應力分布

1.外緣處D點 2.出氣側底部圓角E點 3.內徑處F點(a)壓氣機第十六級輪盤參考點

1.外緣處A點 2.出氣側底部圓角B點 3.內徑處C點(b)透平第一級輪盤參考點圖13 冷態啟動工況轉子上參考點輪盤等效應力隨時間變化曲線

2.2冷態啟動工況的計算結果分析

(1)由于有冷卻空氣的存在,初始階段壓氣機轉子部分的溫度水平較透平的溫度水平高,但隨著負荷的增大,透平轉子的溫度逐漸增大,直至啟動結束,最終透平轉子的最高溫度出現在透平第三級輪盤外緣處,其值達494.97 ℃,而壓氣機轉子的最高溫度出現在壓氣機第十六級輪盤外緣處,其值達402 ℃;溫度基本上沿著轉子的軸向及徑向逐漸增大,其中壓氣機第十七級輪盤由于有轉子冷卻空氣的漏氣冷卻,存在沿軸向負的溫度梯度。

(2)在透平第一級輪盤的出氣側底部圓角處存在很大的應力集中,且在綜合加載應力中離心應力占很大比重,在啟動的整個過程中,最大應力點的分布基本不變,均是在透平輪盤第一級輪盤出氣側底部圓角處,除了幾個應力集中區域的應力較大外,轉子其他部分的綜合加載等效應力大致在0~450 MPa范圍內。

(3)冷態啟動過程中換熱系數的變化按轉速與流量的變化選取,又由于導熱需要時間以及氣流溫度在不同階段的變化,轉子內外溫差在初始階段迅速升高,暖機過程有一定降低,啟動結束階段又有所升高,進入穩定運行階段逐漸減小到一定值;透平轉子中由于有冷卻空氣的存在,整個啟動過程中透平轉子溫度的變化要比壓氣機轉子緩慢。

(4)應力的變化規律與溫度變化規律基本一致,只是應力變化相對溫度變化要有一定的滯后。在輪盤出氣側的底部圓角處,由于應力集中的緣故,應力值均較高。壓氣機轉子的最大綜合加載等效應力的位置是隨啟動過程不斷變化的,參考區域第十六級輪盤的最大等效熱應力出現在輪盤出氣側底部圓角處,并在暖機初始階段內達最大值520 MPa,而在暖機過程中,該處的熱應力有所減小。透平轉子的最大綜合加載等效應力出現在第一級輪盤的出氣側底部圓角處,并在負荷上升完畢后的一段時間內達最大值1116 MPa,該值超過了材料的屈服極限,故該區域是壽命監視的重點區域。

3 熱態啟動工況的計算結果及分析

3.1熱態啟動工況的計算結果

熱態啟動過程中,機組主要經歷了冷加速階段、吹掃階段、點火階段、熱加速階段、自持升速階段及升負荷并網階段。熱態啟動與冷態啟動的主要區別在于啟動的初始時刻轉子的溫度較高,啟動經歷時間短,啟動的溫升率要高于冷態啟動工況,如圖14、圖15所示。

圖14 熱態啟動燃機轉速及功率隨時間變化曲線

1.壓氣機入口 2.壓氣機出口 3.透平出口 4.透平第二級靜葉腔室 5.透平第三級靜葉腔室 6.透平第四級靜葉腔室 7.透平第四級下端腔室 8.轉子冷卻空氣入口圖15 熱態啟動工況燃機各主要監測點溫度隨時間變化曲線

3.2熱態啟動工況的計算結果分析

(1)熱態啟動初始時刻的溫度較高,轉子上最高溫度為190.31 ℃,并且由于各轉子表面氣流溫度的差異較大,雖然換熱系數均較小,但轉子內仍存在較大的溫差;隨著啟動的進行,轉子表面換熱系數及轉速、氣流溫度逐漸增大,轉子溫度場分布逐漸與穩態工況類似,如圖16~圖18所示。

圖16 初始時刻轉子溫度場分布

圖17 啟動后90 min轉子溫度場分布

圖18 啟動后150 min轉子溫度場分布

(2)如圖19所示,熱態啟動參考點溫差變化規律與冷態啟動相似,也有所不同,即:①機組熱態啟動時,轉子金屬的初始溫度較高,氣流溫度低于轉子表面溫度,在壓氣機輪盤上出現了負的溫差,而由于冷卻空氣的存在,透平輪盤上溫度變化平緩,沒有出現負溫差;②透平輪盤內外溫差在暖機過程中并沒有降低,且在暖機結束后的升負荷階段迅速升高到最大值;③整個過程的輪盤內外溫差的變化均沒有冷態啟動過程明顯。

(a)壓氣機第十六級輪盤

(b)透平第一級輪盤1.外緣處D點 2.內徑處F點 3.兩點的溫差圖19 熱態啟動工況輪盤參考點溫度隨時間變化曲線

(3)如圖20~圖22所示,綜合加載等效應力在初始階段以熱應力為主,因而在透平第三級輪盤中部有較大等效應力,而隨著轉速升高,離心應力所占比重上升,其最大應力點仍在透平第一級輪盤出氣側底部圓角處。

圖20 啟動后30 min轉子的綜合加載等效應力分布

圖21 啟動后90 min轉子的綜合加載等效應力分布

圖22 啟動后150 min轉子的綜合加載等效應力分布

(4)如圖23所示,第十六級輪盤的最大綜合加載等效應力出現在輪盤出氣側底部圓角處,并在暖機初始一段時間內達最大值582 MPa,而透平轉子的最大綜合加載等效應力出現在第一級輪盤的中間圓角處,并在升負荷結束后的一段時間達最大值955.30 MPa。

4 兩種啟動工況下的壽命損耗率

通過有限元分析得到參考點的等效應力值,根據Neuber公式及材料的硬化曲線,計算得到該點處的等效真實應變幅為Δε,再由廠方提供的材料低周疲勞壽命曲線計算出危險點的壽命損耗,見表2。

綜上所述,最危險的啟動工況為冷態啟動工況,啟動一次其低周疲勞壽命損耗高達0.0211%。見表3,對照三家主流燃氣輪機制造商的壽命設計數據,F級燃氣輪機的設計標準啟動次數在3200~5000之間[13]。根據GE公司的重型燃氣輪機運行和維護準則,1次冷態啟動可換算為2次標準啟動,1次熱態啟動可換算為0.5次標準啟動。故本文計算的疲勞壽命數據符合設計標準[14]。

1.外緣處D點 2.出氣側底部圓角E點 3.內徑處F點(a)壓氣機第十六級輪盤參考點

1.外緣處A點 2.出氣側底部圓角B點 3.內徑處C點(b)透平第一級輪盤參考點圖23 熱態啟動工況參考點輪盤等效應力隨時間變化曲線

部位等效應變Δε疲勞壽命Nf壽命損耗d(%)冷態啟動壓氣機第十五級輪盤G點0.0028192800.0026壓氣機第十六級輪盤E點0.0023694900.0007透平第一級輪盤B點0.004923640.0211熱態啟動壓氣機第十五級輪盤G點0.0027292150.0017壓氣機第十六級輪盤E點0.0018>1000000<0.0005透平第一級輪盤B點0.003269730.0072

表3 主流燃氣輪機制造商壽命設計標準

5 結論

(1)在本文兩種啟停過程中,冷態啟動是最危險的工況;溫差越大的部位熱應力越大,如透平輪緣處,而綜合加載考慮離心力后最大熱彈性應力值一般都出現在透平第一級輪盤的出氣側底部圓角處及壓氣機后三級輪盤的中間部位的圓角區域,這些都是啟動過程中應該重點關注的地方。

(2)燃機透平轉子中由于有冷卻空氣的存在,能夠降低啟動時因氣流參數迅速變化而導致的轉子內外溫差,其熱應力變化趨于平緩,變化規律與蒸汽輪機轉子有所不同。

(3)轉子應力值的波動取決于轉子表面溫升率及內外表面溫差的變化。暖機過程會使轉子的內外溫差得到緩解,而且對于壓氣機輪盤的緩解較為明顯;相反,氣流參數或機組功率的驟升或驟降會使轉子內外表面溫差快速增大,從而促使轉子的應力值增大。因此,在啟停過程中,操作人員應嚴格控制氣流的參數變化。

(4)基于轉子熱彈性瞬態溫度場及應力場的計算結果對轉子的壽命損耗進行了評估,計算結果表明,冷態啟動工況是最危險的工況,在該工況下透平第一級輪盤的底部圓角處的應力很大,略超過材料屈服極限,其冷態啟動一次的壽命損耗達0.02%,該區域仍是結構形狀優化改進的重點區域,對機組的運行規程進行優化有現實意義。

[1]Antonio C P, Luis S R, Jesus N G,et al.Integration of Thermal Stress and Lifetime Supervision System of Steam Turbine Rotors[C]//American Society of Mechanical Engineers.Berlin,2008: 1035-1044.

[2]Owen J M.Modelling Internal Air Systems in Gas Turbine Engines[J].Journal of Aerospace Power,2006,22(4):515-520.

[3]Owen J M, Wilson M.Some Current Research in Rotating-disc Systems[J].Annals of the New York Academy of Sciences,2001,934(1):206-221.

[4]Owen J M, Rogers R H.Flow and Heat Transfer in Rotating-disc Systems Volume 1-Rotor-Stator Systems[M].Research Studies Press, UK; John Wiley N Y. 1989.

[5]Owen J M, Rogers R H. Flow and Heat Transfer in Rotating-disc Systems Volume 2-Rotating Cavities[M].Research Studies Press, UK; John Wiley N Y. 1995.

[6]曹玉璋,陶智,徐國強,等.航空發動機傳熱學[M].北京:北京航空航天大學出版社,2005.

[7]楊志磊,王海寧,楊承剛,等.600 MW汽輪機轉子熱應力及壽命損耗分析研究[J].汽輪機技術,2011,55(5):383-385.

Yang Zhilei, Wang Haining, Yang Chenggang,et al.Analysis of 600 MW Turbine Rotor Thermal Stress and Loss of Life[J]. Turbine Technology,2011,55(5):383-385.

[8]黎明,楊繼明,白云.利用ALGOR軟件對汽輪機轉子的熱應力分析[J].汽輪機技術,2008,50(4):247-250.

Li Ming,Yang Jiming,Bai Yun.Therma1 Stress Analysis of Steam Turbine Rotor with ALGOR[J]. Turbine Technology,2008,50(4):247-250.

[9]武新華,荊建平,夏松波,等.600 MW汽輪機轉子疲勞壽命計算[J].汽輪機技術,1999,41(3):157-160.

Wu Xinhua,Jing Jianping,Xia Songbo,et al.Fatigue Life Calculation of 600 MW Turbine Rotor[J].Turbine Technology,1999,41(3):157-160.

[10]航空發動機設計手冊總編委會.航空發動機設計手冊[M].北京:航空工業出版社,2001.

[11]張保衡.大容量火電機組壽命管理與調峰運行[M].北京:水利電力出版社,1988.

[12]史進淵,鄧志成,楊宇.超臨界和超超臨界汽輪機轉子葉根槽傳熱系數的計算[J].動力工程學報,2010,30(7):478-484.

Shi Jinyuan, Deng Zhicheng,Yang Yu.Supercritical and Ultra-supercritical Steam Turbine Rotors[J].Journal of Chinese Society of Power Engineering,2010,30(7):478-484.

[13]Swanminathan V P,Dean G J,Scheibel J R.Integreted Approach to Gas Turbine Rotor Condition Assessment and Life Management[C]//Proceedings of ASME Turbo Expo. Copenhagen, 2012:3-4.

[14]Balevic D, Burger R,Forry D.Heavy-duty Gas Turbine Operating and Maintenance Considerations,GER36320K[R].Altlanta:GE Electric Company, 2004.

(編輯陳勇)

Life Estimation for Transient Temperature and Stress Field of Heavy-duty Gas Turbine Rotor

Ou WenhaoYuan QiShi Qingxin

Xi’an Jiaotong University,Xi’an,710049

Taking a specific heavy-duty gas turbine of F class as study object, based on the monitoring data of actual operation, curves denoting temperature changes over time of measuring points in the F class heavy-duty gas turbine unit were extracted at cold starting and hot starting conditions. In these two starting conditions, the heat transfer boundary conditions were calculated throughout, the finite element method was used to analyze the transient temperature field and stress field of the gas turbine rotor in different starting conditions. Based on the results, the low-cycle fatigue life loss was assessed. The calculation results show that maximum thermal elastic stress generally appears at the bottom fillet on the outlet side of the turbine’s first stage disc and rounded area on the middle part of the compressor’s third stage disc; airflow parameter change speed has enormous influences on the rotor’s stress distribution and the life loss of cold starting is up to 0.02%, greater than the corresponding value of the hot starting. The conclusion provides an important reference for the optimization of wheel structure design and the unit operating procedures.

gas turbine rotor; starting condition; finite element method; stress field; life assessment

2014-02-25

TK473DOI:10.3969/j.issn.1004-132X.2015.07.013

歐文豪,男,1989年生。西安交通大學能源與動力學院碩士研究生。主要研究方向為葉輪機械動力學設計及動力學特性。袁奇(通信作者),男,1963年生。西安交通大學能源與動力學院教授、博士研究生導師。石清鑫,男,1987年生。西安交通大學能源與動力學院碩士研究生。

主站蜘蛛池模板: 一级毛片基地| 亚洲三级影院| 无码中文字幕加勒比高清| 最近最新中文字幕在线第一页 | 国产激情在线视频| 国产交换配偶在线视频| 国产在线无码av完整版在线观看| 久久久久免费精品国产| 欧美专区在线观看| 国产91熟女高潮一区二区| 国产在线视频欧美亚综合| 日本尹人综合香蕉在线观看| 99re免费视频| 免费人欧美成又黄又爽的视频| 国产视频久久久久| 天天躁狠狠躁| 一级毛片在线播放免费| 亚洲AV色香蕉一区二区| 国产精品久久久久久久久久久久| 国产电话自拍伊人| 亚洲国产成熟视频在线多多| 老司机久久99久久精品播放| 综1合AV在线播放| 亚洲国产精品日韩欧美一区| 成人av手机在线观看| 麻豆AV网站免费进入| 亚洲资源站av无码网址| 国产97色在线| 99无码中文字幕视频| 国产成人AV大片大片在线播放 | 国产成人91精品| 2020国产在线视精品在| 精品国产成人三级在线观看| 欧美亚洲国产精品第一页| 国产精品亚洲а∨天堂免下载| 影音先锋亚洲无码| 亚洲an第二区国产精品| 国产色图在线观看| 人人看人人鲁狠狠高清| 中国一级毛片免费观看| 国产精品55夜色66夜色| 18禁黄无遮挡免费动漫网站| 国产中文在线亚洲精品官网| 国产97公开成人免费视频| 日日噜噜夜夜狠狠视频| 老司机久久99久久精品播放| 欧美日韩午夜| 91精品啪在线观看国产| 亚洲欧美日韩另类在线一| 98精品全国免费观看视频| 国产h视频在线观看视频| 理论片一区| 尤物精品视频一区二区三区| 久久久久久久蜜桃| 99精品视频在线观看免费播放| 国产亚洲精品91| 超薄丝袜足j国产在线视频| 成人午夜视频网站| 色综合天天视频在线观看| 日韩a级片视频| yjizz国产在线视频网| 亚洲三级a| 欧美性猛交一区二区三区| 国产色网站| www.91中文字幕| 米奇精品一区二区三区| 久久久久亚洲精品成人网| 综合久久久久久久综合网| 人妻少妇乱子伦精品无码专区毛片| 99视频只有精品| 国产精品视频观看裸模| 亚洲国产综合自在线另类| 国产精品嫩草影院视频| 国产精品免费电影| 在线观看无码av免费不卡网站| 91久久青青草原精品国产| 亚国产欧美在线人成| 无码AV日韩一二三区| 亚洲婷婷丁香| 国产成人综合在线视频| 欧美在线观看不卡| 国产尤物jk自慰制服喷水|