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

變工況條件下汽輪機高壓缸末級氣動及強度性能研究

2023-01-10 06:51:30石紅暉王海波曹蓉秀姚力晏鑫
發電技術 2022年6期
關鍵詞:模型

石紅暉,王海波,曹蓉秀,姚力,晏鑫*

(1.國家能源集團科學技術研究院有限公司,江蘇省 南京市 210046;2.西安交通大學能源與動力工程學院,陜西省 西安市 710049)

0 引言

對大中型火電機組進行供熱改造,可以在不增加電力裝機規模的前提下,發展熱電聯產集中供熱,該措施既能解決分散小鍋爐供熱帶來的能源、環境問題,又能滿足日益增長的居民采暖和工業用汽需求[1-2]。然而,抽汽供熱時主蒸汽流量的變化將對汽輪機高壓缸末級的運行狀態、氣動性能以及葉片、輪盤等結構的強度造成一定影響[2]。因此,研究不同抽汽量工況下高壓缸末級的氣動和強度性能,對于汽輪機供熱改造方案的制定以及熱電聯產機組在變工況條件下的安全、穩定、高效運行具有重要的指導意義。

目前,關于汽輪機供熱改造變工況特性研究主要集中在小流量工況下低壓缸的流動特性[3-5],以及低壓缸末級長葉片的強度和振動方面[6-8]。然而,與低壓缸末級長葉片不同的是,由于高壓缸葉片相對較短,葉片所受離心力和氣動力的相對作用強度發生明顯變化;且短葉片多采用T型葉根,其受力特點、結構強度與長葉片普遍采用的樅樹型葉根差別較大[9-10];另外,固定在動葉頂部的汽封結構也承受一部分氣動力,而低壓缸葉片的強度分析一般較少考慮汽封的影響。因此,研究帶T型葉根、葉頂汽封的高壓缸葉片在離心力、氣動力作用下的應力狀態和變形情況十分必要。

對于帶T型葉根葉片的強度分析,部分學者進行了研究。Shankar等[11]建立了葉片-輪盤三維有限元模型,采用葉根與輪盤接觸面節點完全匹配的六面體網格以提高計算精度,研究結果表明,在離心力的作用下,T型葉根進汽側下倒圓處應力最大,葉型壓力面根部的應力次之;其不足之處在于將葉根與輪盤視作完全緊密配合,忽略了二者之間非接觸面的安裝間隙。Gowda等[12]考慮了葉根與輪盤的實際配合間隙,采用旋轉周期對稱邊界研究了整圈葉片的強度和疲勞特性,結果表明,葉根徑向接觸面圓角處等效應力最大且易出現疲勞損傷。文獻[13-15]在T型、雙T型葉根和葉根槽的結構設計和優化方面開展了較為充分的研究,對葉根和輪盤接觸過程的有限元分析具有一定的指導作用。另外,Ehlers等[15]提出的精確適應復雜幾何結構的多塊六面體網格劃分方法,對于提高有限元分析精度、節省計算資源都十分有利。然而,上述研究主要考慮離心力對葉片強度的影響,對葉片在三維流場中的受力情況以及變工況運行時的強度安全特性未進行深入研究。

為了研究變工況條件下高壓缸末級的氣動和強度性能,本文建立了帶汽封結構的高壓缸末兩級計算流體動力學(computational fluid dynamics,CFD)計算模型,以及帶外包T型葉根和葉頂汽封的末級動葉-輪盤的有限元分析(finite element analysis,FEA)模型,采用單向流固耦合[16]分析方法,對比了不同抽汽量(0%、8%、15%和20%)工況下高壓缸末兩級的總體性能和主要氣動參數沿葉高的分布,分析了離心力和氣動力對葉片?輪盤系統靜應力和變形位移的影響機制,揭示了葉片、輪盤關鍵部位應力和葉頂最大變形位移隨抽汽量的變化規律。

1 數值計算方法

1.1 CFD計算模型及方法

圖1為帶汽封結構的高壓缸末兩級(第11、12級)CFD計算模型,表1為葉柵和汽封的主要幾何參數。由于兩級動、靜葉片數分別為64和80,故計算模型設置為4個動葉通道和5個靜葉通道相匹配,使得動、靜計算域周向寬度相同,從而其信息傳遞可以采用凍結轉子法(frozen rotor),以減小數值計算誤差[17]。汽封與葉柵通道周向寬度保持一致,均為22.5°。

圖1 高壓缸末兩級CFD計算模型Fig.1 CFD computational model for the last two stages in high pressure cylinder

表1 葉柵和汽封的主要幾何參數Tab.1 Main geometrical parameters for blades and seals

圖2為葉柵和汽封的網格模型。利用商業軟件NUMECA-AutoGrid5生成兩級葉柵的結構化網格,網格最小偏斜角為28.9°;利用軟件ANSYSICEM生成汽封的結構化網格,網格最小偏斜角為40.5°。葉柵和汽封近壁面第1層網格高度均為0.001 mm,網格節點總數約為5 073萬。

圖2 葉柵和汽封的計算網格Fig.2 Computational meshes for blades and seals

CFD計算模型的子午面視圖如圖3所示,其中:SD表示靜止域;RD表示旋轉域,其轉速為3 000 r/min。表2為設計工況下CFD計算模型的邊界條件。模型進口給定總溫、總壓,軸向進汽,中等湍流強度(5%),出口給定平均靜壓。壁面采用絕熱、光滑、無滑移邊界,汽封中與輪轂面、動葉相連的壁面為旋轉邊界,轉速為3 000 r/min;與圍帶面、靜葉相連的壁面為靜止邊界;模型周向2個側面采用旋轉周期邊界。動、靜葉柵之間,以及汽封與主通道之間交界面的位置、編號如圖3所示。動、靜計算域(編號(1)—(11)所示)之間信息傳遞采用凍結轉子法[18]。

圖3 CFD計算模型子午面視圖Fig.3 Meridian view of CFD computational model

表2 CFD計算模型的邊界條件Tab.2 Boundary conditions of CFD computational model

利用軟件ANSYS-CFX求解三維定常RANS方程,模型工質為IAPWS-97標準的可壓縮水蒸汽,湍流模型選用SST k-ω兩方程湍流模型,控制方程離散采用高分辨率格式。當控制方程的均方根殘差小于1×10?5,相鄰迭代步出口截面平均溫度波動小于0.01 K時,認為計算收斂。

為得到網格無關解,建立了5套近壁面第1層網格高度ywall的葉柵和汽封網格,采用動、靜葉柵通道數為1∶1的單流道計算模型,研究了ywall對質量流量m和透平級總-總等熵效率ηit的影響,計算結果如表3所示。其中,透平級總-總等熵效率定義為

表3 近壁面第1層網格距離無關性分析Tab.3 Independency analysis of the first mesh layer distance near wall

式中:i0*為透平級進口總焓;i2*為透平級出口總焓,i2's為透平級出口等熵靜焓。

由表3可知,兩級總?總等熵效率對ywall不敏感,而質量流量隨ywall的減小而減小。當ywall由0.001 mm減小到0.000 5 mm時,流量變化已不再明顯,故最終采用ywall=0.001 mm的網格,此時近壁面平均y+=7.8,滿足SST k-ω兩方程湍流模型的計算要求。

針對葉柵和汽封,研究了網格密度對m和ηit的影響,計算結果分別如表4、5所示。可見,網格密度對計算結果的影響不顯著。最終確定葉柵和汽封網格節點數分別為710萬和479萬,即單流道網格總節點數為1 189萬。

表4 葉柵網格密度無關性分析Tab.4 Independency analysis for grid density of blades

表5 汽封網格密度無關性分析Tab.5 Independency analysis for grid density of seals

1.2 FEA計算模型及方法

帶外包T型葉根和葉頂汽封的高壓缸末級動葉-輪盤FEA計算模型如圖4所示。葉型高度為118.3 mm,葉型根部輪轂面半徑為450 mm。取整周葉輪的1/64扇區進行單葉片的強度分析,模型周向寬度為5.625°。葉片和輪盤材料均選用2Cr13,其在400℃下的主要力學性能參數[19]如表6所示。由CFD計算結果可知,設計工況下高壓缸末級動葉工作溫度約為400℃,變工況時溫度最高上升了約10℃。由于溫度變化對材料的力學性能影響微小,故采用400℃下材料的力學性能進行計算,忽略溫度變化引起的熱應力對葉片強度的影響。

圖4 高壓缸末級動葉-輪盤FEA計算模型Fig.4 FEA computational model for the last stage rotor blade-disk in high pressure cylinder

表6 2Cr13在400℃下的主要力學性能參數Tab.6 Main mechanical property parameters of 2Cr13 at 400℃

圖5葉片?輪盤FEA計算網格Fig.5 FEA computational meshes for blade-disk

圖5 為葉片-輪盤FEA網格模型。為了提高計算精度、節省計算資源,采用分塊劃分網格的方法,除輪盤軸線處存在少量楔形單元外,模型其余部位均為六面體單元;葉根與輪盤接觸面(標號A、B、C、D),以及葉根、葉頂汽封與葉型連接面(標號E、F)網格節點一一對應;對接觸面附近和其他高應力區網格進行加密,最小單元尺寸為0.3 mm。由于CFD網格與FEA網格節點不完全對應,為了減小氣動力耦合時帶來的插值誤差[16],將氣動力耦合面網格尺寸設置得盡可能小,最大單元尺寸為2 mm。計算單元采用適用于接觸分析問題的C3D8I(8節點非協調線性六面體)單元,其總數為34.8萬。

圖6為FEA計算模型的邊界條件及載荷。對輪盤旋轉軸所在邊線設置位移全約束,忽略主軸沿軸向的滑動,對輪盤主軸部分軸向2個側面設置軸向位移約束。由于整圈葉片沿輪周方向安裝,葉根、葉頂汽封處于互相擠壓的緊密接觸狀態,故在葉根、葉頂汽封和輪盤的周向2個側面設置切向位移約束。需要說明的是,由于單個扇區的氣動力載荷沿周向分布不具有周期性,因此計算模型周向2個側面無法采用旋轉周期對稱邊界,而切向位移約束是相對合理的[12]。為了較為準確地模擬葉根與輪盤槽之間的相互作用,將軸向接觸面A、B和徑向接觸面C、D設置為面面接觸;接觸特性設置為小滑移;接觸面法向行為設置為硬接觸,切向行為忽略摩擦作用[8]。對整個模型施加離心力載荷(如圖6中綠色箭頭所示),角速度ω=100π rad/s。為了準確模擬葉片在三維流場中的受力情況,將CFD計算得到的葉型表面、輪轂面、圍帶面和葉頂汽封表面的氣動力耦合到FEA計算模型上,4個部位的氣動力載荷如圖6中紫色箭頭所示。

圖6 FEA計算模型邊界條件及載荷Fig.6 Boundary conditions and loads for the FEA computational model

2 計算結果分析

2.1 變工況透平級總體參數

表7為設計工況(0%抽汽)及抽汽量占設計工況總流量的8%、15%、20%工況下高壓缸末兩級總體參數。抽汽口位于高壓缸次末級(第11級)進口處,忽略周向壓力分布不均的影響。次末級進口總壓Pin作為邊界條件給定,隨抽汽量的增大近似呈線性減小。各工況的進口總溫和出口靜壓保持一致,分別為708.05 K和7.76 MPa。隨著抽汽量的增大,出口總溫Tout逐漸升高,20%抽汽量工況下Tout較設計工況上升了約10 K,高壓缸末級葉片工作環境變差;兩級輸出功率P近似呈線性減小,20%抽汽量工況下P較設計工況降低了約44%。抽汽供熱使高壓缸末兩級對外做功減少,從而影響機組運行的經濟性。

表7 高壓缸末兩級變工況總體參數Tab.7 Overall parameters of final two stage variable conditions of high-pressure cylinder

2.2 變工況透平級主要氣動參數

為了說明汽封泄漏引起的二次流對透平級氣動參數的影響,圖7給出了20%抽汽量工況下葉頂和葉根附近二次流的子午面流線。由圖7(a)可知,汽封R11的泄漏流進入末級葉柵通道后對主流產生擾動并不斷向后擴展,二次流影響范圍為葉頂附近約20%葉高區域;汽封R12對末級流動影響較小。由圖7(b)可知,汽封S11的泄漏流對后三排葉柵的流動均產生影響,發展到末級動葉時二次流影響范圍為葉根附近約30%葉高區域;汽封S12與S11的泄漏流在末級動葉相互摻混,其影響范圍為葉根附近約20%葉高區域。

圖7 20%抽汽工況下汽封泄漏引起的二次流子午面流線Fig.7 Secondary flow meridian streamlines caused by steam steal leakage under 20%extraction condition

圖8給出了不同抽汽量工況下兩級總-總等熵效率ηit隨相對葉高Hˉ的變化規律。對于第11級,ηit在30%葉高以上區域變化較為平緩,且隨抽汽量的增大而減小,20%抽汽量工況下ηit較設計工況減小了約3%;30%葉高以下區域由于受汽封S11泄漏流的影響,曲線在20%葉高處出現極值點,透平級氣動效率最差,且該極值點隨抽汽量的增大逐漸向右下方移動,其原因是汽封泄漏量隨主流量的減小而減小,使得二次流的影響區域縮小,效率最低點向葉根方向移動。對于第12級,在30%葉高以上區域,抽汽量較小時ηit與設計工況較為接近,20%抽汽量工況下ηit偏離設計工況較多;30%葉高以下區域由于同時受汽封S11和S12泄漏流的影響,曲線發生明顯波動,且變工況時ηit嚴重偏離設計工況,10%葉高處最大偏離量約為3.8%。可見,抽汽供熱對高壓缸末兩級的氣動效率影響顯著。

圖8總-總等熵效率沿葉高的分布Fig.8 Total-total isentropic efficiency distributions along spanwise direction

圖9 反動度沿葉高的分布Fig.9 Reaction degree distributions along the spanwise direction

式中:i1為動葉進口靜焓;i2s為動葉出口等熵靜焓。

由圖9可知,在30%葉高以上區域,兩級反動度Ω沿葉高方向逐漸增大,由于受汽封泄漏流的影響,葉根附近Ω出現波動。在第11級動葉50%葉高以下區域,Ω隨抽汽量的增大而增大,20%抽汽量工況下Ω較設計工況最大增加了約0.9%;葉頂附近Ω幾乎不受抽汽影響。第12級Ω沿葉高的分布呈現出類似的規律,但變工況時Ω偏離設計工況更為明顯,20%抽汽量工況下Ω較設計工況最大增加了約1.6%,且該極值點也出現在20%葉高處。可見,抽汽供熱對高壓缸末級50%葉高以下區域的反動度影響較大,使得末級動、靜葉柵焓降分配發生變化,從而影響機組運行的經濟性和安全性[20]。

圖10給出了不同抽汽量工況下兩級葉柵出口汽流角隨相對葉高Hˉ的變化規律。其中,靜、動葉柵出口汽流角分別定義如下:

式中:Vz為絕對汽流速度軸向分量;Vt為絕對汽流速度切向分量;Wz為相對汽流速度軸向分量;Wt為相對汽流速度切向分量。

從圖10可以看出,四排葉柵葉展中部汽流角變化較平緩,而靠近葉頂和葉根區域由于受汽封泄漏流的影響,汽流角發生明顯波動,尤其是葉根附近。動、靜葉柵30%~80%葉高區域汽流角基本不隨抽汽量變化。除第11級靜葉外,變工況時其余3排葉柵的葉頂和葉根附近出口汽流角明顯偏離設計工況,且抽汽量越大,偏離越嚴重,尤其是20%和90%葉高處。20%抽汽量工況下第12級動葉20%葉高出口汽流角較設計工況最大增加了約2.4°。抽汽供熱對高壓缸末級葉柵出口汽流 角的影響大于前一級,且對動葉的影響大于靜葉。葉柵出口汽流角的分布規律基本與圖7所示二次流的影響范圍相對應。與設計工況相比,抽汽工況下葉柵出口汽流角的變化容易引發較強的流動分離和尾跡渦,使流動損失增大、內效率降低[18]。

圖10 葉柵出口汽流角沿葉高的分布Fig.10 Blade outlet flow angle distributions along the spanwise direction

2.3 離心力和氣動力對葉片強度的影響對比

圖11為離心力單獨作用時葉片和輪盤的應力云圖。由于整個葉片重心與輪盤中心線不完全重合,輪盤左右兩側應力分布并非完全對稱,最大應力點位于葉根槽進汽側倒圓(輪盤與葉根徑向接觸面C的根部倒圓)處,應力為428.9 MPa。相應地,葉片最大應力點位于葉根進汽側下倒圓(葉根與輪盤徑向接觸面C的根部倒圓)處,該結果與文獻[11-12]結論一致。葉根進汽側下倒圓處的最大應力為423.3 MPa,約為出汽側下倒圓處最大應力的1.3倍。葉型表面最大應力點位于壓力面根部40%軸向弦長處,應力為286.6 MPa,而吸力面應力相對較小。由于材料屈服極限為735 MPa,故在離心力作用下葉片和輪盤強度滿足要求,安全系數為1.74。

圖11 離心力單獨作用時葉片和輪盤應力分布Fig.11 Stress distributions in blade and disk with only centrifugal force

圖12給出了設計工況下FEA計算模型中氣動力耦合面的壓力分布與CFD計算結果的對比。由圖可知,FEA壓力云圖的色帶邊緣為鋸齒狀,這是氣動力耦合面在數據傳遞中進行插值運算的結果[16],可通過加密FEA模型中氣動力耦合面附近的計算網格,以減小插值誤差。總體來看,FEA計算模型中4個氣動力耦合面的壓力分布與CFD計算結果吻合良好。

圖12 設計工況下FEA和CFD交界面上參數分布Fig.12 Parameter distributions on the interface between FEA and CFD domain under design condition

圖13為設計工況下離心力、氣動力作用時葉片和輪盤的應力云圖。通過對比圖11、13可知,氣動力使葉片和輪盤進汽側應力遠遠高于出汽側,這是由汽流沿軸向的壓力差所決定的。葉片最大應力點依然位于葉根進汽側下倒圓處(該結果與文獻[14]結論相一致),應力為420.7 MPa,葉根槽進汽側倒圓處的應力為379.5 MPa,二者較離心力單獨作用時分別降低了約0.6%和11.5%。葉根進汽側下倒圓處的最大應力約為出汽側下倒圓處最大應力的3.4倍,而這一比例在離心力單獨作用時僅為1.3;葉型壓力面根部最大應力點位置基本不變,應力由286.6 MPa增大為289.1 MPa。由此可見,離心力對葉片主要產生沿徑向的拉伸作用,而氣動力對葉片主要產生沿軸向的彎曲作用,使葉片受力更不均勻;氣動力的切向彎曲作用相對較弱。

圖13 離心力、氣動力作用時葉片和輪盤應力分布Fig.13 Stress distributions in blade and disk with both centrifugal force and aerodynamic force

圖14為加載氣動力前后葉片變形位移云圖(放大50倍)。無氣動力作用時,葉頂最大位移為0.216 mm,主要為徑向位移。加載氣動力后,葉頂最大位移增大為原來的2倍以上,為0.443 mm,徑向位移基本不變,而軸向位移明顯增大。輪盤和葉根進汽側位移明顯大于出汽側,葉型和葉頂汽封位移分布沿徑向和軸向的變化速度均明顯增大。由此可見,氣動力對葉片變形位移影響顯著,故在葉片強度分析中需同時考慮離心力和氣動力的作用,且不能只關注最大應力的變化。

圖14 加載氣動力前后的葉片變形位移分布(放大50倍)Fig.14 Deformation displacement distributions of blade before and after loading aerodynamic force(magnify 50 times)

2.4 變工況末級動葉強度分析

圖15 為葉片和輪盤關鍵部位高應力區的位置及編號,表8為不同抽汽量工況下相應位置的應力。可以看出,在離心力和氣動力的作用下,葉根最大應力點位于進汽側下倒圓處(位置⑤);葉型表面最大應力點位于壓力面根部40%軸向弦長處(位置②);輪盤最大應力點位于葉根槽進汽側倒圓處(位置⑥),在0%、20%抽汽量工況下,位置⑥處應力較位置⑤處應力分別降低了10%、6%。隨著抽汽量的增大,除了葉頂汽封與葉型吸力面連接處(位置①)的應力略微增大(不超過設計工況應力值的0.2%)之外,其余部位應力均近似呈線性減小;20%抽汽量工況下葉根最大應力較設計工況減小了約7.1%。

圖15 葉片、輪盤高應力區位置及編號Fig.15 Position and number of high stress area in blade and disk

表8 變工況葉片、輪盤關鍵部位的應力Tab.8 Stress of key parts of blade and disk under variable conditions MPa

表9給出了不同抽汽量工況下葉片頂部最大變形位移。結合圖14(b)可知,葉片變形位移沿葉高向上近似呈線性增大,且進汽側位移明顯大于出汽側。葉片變形主要沿軸向和徑向,徑向位移主要由離心力的拉伸作用而產生,由于轉速保持不變,變工況時葉片徑向位移變化較小(不超過設計工況徑向位移的0.3%);軸向位移主要由氣動力的彎曲作用而產生,由于末級動葉前后壓力差隨抽汽量的增大而減小,氣動力的彎曲作用減弱,軸向位移隨之減小。葉頂最大變形量隨抽汽量的增大近似呈線性減小,20%抽汽量工況下最大變形量較設計工況減小了約19.2%。

表9 變工況葉片頂部最大變形位移Tab.9 Maximum deformation displacement of blade tip under variable conditions mm

3 結論

采用單向流固耦合分析方法,對4種抽汽量工況下汽輪機高壓缸末級的氣動和強度性能進行了數值研究,主要結論如下:

1)隨著抽汽量的增大,高壓缸出口總溫逐漸升高,應注意末級葉片的運行安全;高壓缸總?總等熵效率顯著下降,輸出功率隨著抽汽量的增大而近似呈線性減小。

2)隨著抽汽量的增大,葉片根部反動度逐漸增大,但葉頂反動度幾乎不受影響;葉展中部汽流角變化較小,但葉頂和葉根區域由于受汽封泄漏流的影響而變化顯著。

3)隨著抽汽量的增大,除葉頂汽封與葉型吸力面連接處的應力略微增大外,其余部位的應力近似呈線性減小;葉根最大應力和葉頂最大變形量均隨抽汽量的增大近似呈線性減小。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 玖玖精品视频在线观看| 亚洲国产综合自在线另类| 视频一区视频二区日韩专区| 国产不卡一级毛片视频| 东京热av无码电影一区二区| 久久性视频| 亚洲制服丝袜第一页| 国产美女主播一级成人毛片| 国产原创演绎剧情有字幕的| 国产95在线 | 美女一区二区在线观看| 中国毛片网| 五月天久久婷婷| 久久99国产综合精品女同| 四虎免费视频网站| 亚洲美女操| 精品国产黑色丝袜高跟鞋| 日韩在线1| 国产乱人伦精品一区二区| 国产sm重味一区二区三区| 91一级片| 亚洲天堂网视频| 亚洲成a人片7777| 久久精品中文字幕少妇| 免费观看成人久久网免费观看| 亚洲最黄视频| 国产黄色片在线看| 国产一级裸网站| 九九久久99精品| 一本大道东京热无码av| 成人午夜视频免费看欧美| 亚洲黄色成人| 69av在线| 狼友视频一区二区三区| 尤物成AV人片在线观看| 国产主播一区二区三区| 精品亚洲国产成人AV| 中文字幕永久在线看| 91国内在线观看| 国产精品亚洲精品爽爽| 色综合成人| 在线精品亚洲国产| 欧美人与性动交a欧美精品| 精品人妻一区无码视频| 漂亮人妻被中出中文字幕久久| 亚洲综合18p| 欧美精品一二三区| 亚洲无码视频一区二区三区| 这里只有精品在线播放| 青青操国产视频| 国产传媒一区二区三区四区五区| 亚洲人成网7777777国产| 制服丝袜亚洲| 在线无码九区| 亚洲日韩国产精品综合在线观看| 国产精品专区第1页| 欧美区一区| 亚洲精选无码久久久| 欧美在线精品怡红院| 亚洲天堂区| 在线免费看片a| 亚洲欧洲日产国产无码AV| 国产在线视频欧美亚综合| 国产91在线免费视频| 久久国产精品电影| 日韩精品一区二区三区免费在线观看| 伊人婷婷色香五月综合缴缴情| 精品少妇人妻av无码久久| 97国产精品视频自在拍| 天天激情综合| 18黑白丝水手服自慰喷水网站| 特级欧美视频aaaaaa| 国产主播一区二区三区| 免费日韩在线视频| 国产福利影院在线观看| 国产主播一区二区三区| 国产在线视频福利资源站| 成人毛片在线播放| 秋霞国产在线| 人妻免费无码不卡视频| 精品夜恋影院亚洲欧洲| 欧美伦理一区|