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

自然循環下鈉冷快堆中迷宮型節流件水力特性研究

2023-12-16 05:20:50周志偉陸道綱
核科學與工程 2023年5期
關鍵詞:實驗模型

周志偉,曹 鋒,陸道綱,曹 瓊

自然循環下鈉冷快堆中迷宮型節流件水力特性研究

周志偉1,曹鋒2,陸道綱2,曹瓊2,*

(1. 中國原子能科學研究院,北京 102413;2. 華北電力大學 非能動核能安全技術重點實驗室,北京 102206)

迷宮型節流件常被安裝在鈉冷池式快堆主容器冷卻系統、泵冷卻系統等流動旁路中,使得進入這些旁路的流量達到設計值,來有效冷卻設備,保證快堆的安全。目前國內對自然循環條件下節流件的阻力特性研究較少,本文對迷宮節流件在自然循環和強制循環兩種條件下進行了實驗和數值研究。首先,實驗采用恒定溫度(84 ℃)的液態水作為工質,得到了兩種循環條件下節流件的阻力系數;其次,對節流件進行數值模擬,將不同湍流模型、壁面法向網格間距下的數值模擬結果與實驗結果進行了比較,驗證了數值模擬的有效性;最后,通過比較不同湍流模型的計算結果和流場細節,得出適合于兩種條件下迷宮式節流件數值模擬的湍流模型,并發現回流漩渦是節流件的主要耗能方式,兩種條件下模擬結果的差異主要集中在二次回流區域大小的不同,自然循環下回流區域明顯大于強制循環,且收縮速度較慢。

迷宮型節流件;水力特性;湍流模型;數值模擬

鈉冷快堆(SFR)具有高效利用鈾資源、嬗變長壽命錒系元素和固有安全性等突出優勢,是第四代先進核反應堆堆芯之一[1]。在SFR中,迷宮式節流件常常安裝在SFR主容器冷卻系統、泵冷卻系統等流動旁路中,以將旁路中的流量和壓降控制在設計范圍內,防止出現旁路中流量過大而導致主回路中冷卻劑流量不足的事故。因此,準確預測節流件的壓降和水力特性對鈉冷快堆的安全評估具有重要意義。

Rhode和Sobolik使用有限差分法模擬了矩形齒迷宮密封中可壓縮流體的泄漏流[3]。Yang和Luis量化了進口圓角對水潤滑密封件的影響。發現彎曲進口密封的進口壓降較小,但由于固體顆粒的侵蝕,進口邊緣容易磨損[4]。Artur等人優化了帶有兩個斜翼的直通式迷宮節流件。同時將數值計算結果與實驗數據進行了比較,得出位置、傾角和翅片厚度對減少泄漏的影響最大[5]。Zhang等人討論了渦流制動器對迷宮密封流場和轉子動力學性能的影響。結果表明,渦流制動器數量、長度的增加,以及渦流制動器間隙的減小,都會改善直線渦流制動器的轉子動力學性能[6]。Wang等人研究了高壓控制閥盤中由一系列直角旋轉組成的曲折迷宮通道。研究發現,與“平行通道”相比,“串聯通道”可以更有效地降低壓降[7]。Asok等人分別研究了圓形腔、方形腔和彎曲腔迷宮密封的流動特性。發現空腔內的多個再循環區域促進了動量耗散,從而提高了密封性能[8]。Wang等人分析了軸向位置對階梯迷宮密封密封性能的影響。結果表明,泄漏和傳熱特性實際上與軸向位置無關[9]。秦亥琪等對迷宮密封進行了實驗和數值研究,得到了其水力特性隨密封幾何尺寸的變化規律[10]。宋怡等人對迷宮式節流件進行了數值分析,并擬合了壓降與幾何參數之間的經驗關系式[11]。除了數值模擬,關于節流密封的實驗也有很多。Dobrowolski等人采用微粒圖像測速儀(PIV)觀察迷宮密封中的潤滑脂流動,評估了迷宮密封中的潤滑脂速度分布,發現潤滑脂稠度與旋轉環轉移速度之間存在關系[12]。Untarou等人進行了數值和實驗研究,以計算固定直通迷宮密封的流動參數。同時清楚地觀察到迷宮密封中的局部流動現象[13]。

綜上所述,國內外有許多針對迷宮密封這類細微通道內流動特性的研究,但對節流件這種安裝在旁路中的節流元件研究較少。在對反應堆進行熱工水力數值仿真時,研究人員往往更關心反應堆堆芯及主回路管道內的流動及傳熱特性,對節流件這類旁路中的節流元件往往進行忽略或簡化。但在自然循環條件下,其阻力特性對整個回路的影響不能忽視。

因此本文將通過實驗和數值模擬來研究迷宮節流件在自然循環和強制循環條件下的阻力特性。首先根據實驗得到其阻力系數。然后采用不同的湍流模型對節流件進行數值模擬,通過與實驗結果進行比較,得出合適的湍流模型。最后,分析節流件在兩種條件下流場的差別。

1 迷宮式節流件的水力實驗

1.1 迷宮節流件結構及尺寸

三維迷宮式節流件由厚度相同、高度相同且高于管道半徑的節流板組成。相鄰節流板之間的間距相同,并以180°交錯排列。節流部分置于圓管內,示意圖如圖1所示。

圖1 迷宮型節流件示意圖[11]

本文研究的節流件由七個節流板組成,兩個相鄰的節流板是交錯的。其具體幾何參數如表1所示。橫截面如圖2所示。

表1 節流件幾何參數

圖2 節流件橫截面

1.2 實驗原理及參數

本實驗是研究迷宮節流件流動特性的驗證實驗。考慮到實驗的安全性,本實驗用去離子水代替液態鈉。為了保持SFR的水力特性,必須完成液態鈉和水之間實驗參數的相似轉換。水力相似轉換需要滿足幾何相似性、運動相似性和動力相似性。

為了滿足幾何相似性,實驗的節流件與鈉冷快堆中的迷宮節流件具有完全相同的幾何參數。同時,由于液態鈉和去離子水都被認為是不可壓縮的粘性流體,因此粘性力是控制其流場水力特性的主要作用力。當兩個流場的雷諾數相同時,保證了運動相似性和動力相似性。雷諾數()是一個無量綱數,表示粘性流體的慣性力與粘性力之比。

式中:——雷諾數的特征長度,在實際計算中選用節流件的內徑;

——流體特征速度,計算中選取入口流速;

鈉冷快堆中流經迷宮節流件的液態鈉溫度為360 ℃。鈉的物性參數計算公式如表2所示。由于實驗的節流件尺寸完全相等,且在相同的單值條件下,要滿足雷諾數相等,即

表2 液態鈉物性[14]

在迷宮節流件中,節流板引起的局部阻力是節流件壓降的主要來源,根據局部阻力系數公式,管道中流動的壓降與阻力系數滿足如下關系式:

為局部阻力系數,是雷諾數的單值函數,即

因此,當不同流體的雷諾數相同時,它們的阻力系數也相同。當獲得去離子水流經被測節流件產生的壓降時,便可根據歐拉相似準則計算液態鈉流經同一節流件產生的壓降。

歐拉相似準則:

因此,本實驗僅需測量被測節流件兩側的水流量和壓降。相似轉換后的迷宮式節流件水力參數如表3所示。

表3 節流件水力參數

1.3 實驗裝置及流程

整個實驗系統由實驗回路和加熱回路兩部分組成。實驗回路包括去離子水箱及其出口閥、被測節流件、流量計、壓力表、節流進口電動閥、頂部排氣閥、變頻泵和熱電偶。加熱回路包括一個額定功率為120 kW(380 V)的筒式加熱器,通過入口和出口閥與實驗回路相連。系統圖如圖3所示。

圖3 實驗回路示意圖

試驗開始時,打開兩個回路中的閥門,變頻泵將去離子水從去離子水箱泵入試驗系統,直到去離子水完全充滿兩個回路,然后打開試驗回路頂部的排氣閥,排出多余的空氣。然后,通過水箱出口閥將去離子水箱與測試系統隔離,以防止實驗過程中加熱的水返回水箱。最后,變頻泵驅動去離子水連續通過加熱回路和實驗回路。由于此時加熱回路與實驗回路相連,整個實驗系統中的水由圓柱形加熱器加熱,直到達到實驗溫度要求。去離子水達到所需溫度后,關閉加熱器,并通過關閉加熱回路的入口和出口閥將實驗回路與加熱回路隔離。整個實驗回路全部外覆保溫材料,保證在實驗過程中去離子水溫度保持恒定。

通過改變電動閥的開度,可以精確控制去離子水的流量。被測節流件的額定流量為11.5 kg/s,試驗期間進行了兩種工況:一種是額定流量的2.5%~6.2%的自然循環工況,另一種是額定流量的75%~110%的強制循環工況。

鈉冷快堆正常運行時,通過主泵帶動鈉冷卻劑流動,為強制循環工況。在鈉冷快堆停堆時,鈉冷卻劑依然存在正常流量6%左右的流量來導出堆芯的余熱,這種情況即為鈉冷快堆的自然循環工況,本實驗的試驗臺架則通過控制去離子水的流量來模擬鈉冷快堆的這兩種循環工況。

當流量穩定時,由計算機采集流量和壓降數據。每個工況持續3~4 min,測量三次。最終結果取平均值。

1.4 實驗結果和不確定度

實驗測量了5種自然循環工況和12種強制循環工況下被測節流件的流量和壓降。結果如圖4所示。

圖4 壓降隨質量流量變化的實驗結果

實驗不確定性有兩個來源,包括數據處理和測量設備。對于鈉冷卻劑,壓降和流量的不確定度分別小于1%和3%。熱電偶已校準,最大不確定度為±0.5 K。

2 迷宮節流件的數值模擬

2.1 模型和網格化方法

數值模擬的幾何模型的大小與實驗1:1相同。采用ICEM軟件劃分非結構化四面體網格,在節流板壁面生成5層邊界層網格,靠近壁面的第一個邊界層網格厚度為9.7×10-5m。ICEM中節流板處的網格如圖5所示。

圖5 節流件網格示意圖

2.2 邊界條件和湍流模型

表4展示了數值模擬的邊界條件。自然循環和強制循環條件下水力試驗和數值模擬的范圍分別為8 124~19 889和244 323~352 796。因此,模擬的液態鈉流動屬于高雷諾數條件下的湍流。采用雷諾平均(RANS)方法下的五種湍流模型。包括標準-、Realizable-、標準-、Spalart-Allmaras(S-A)和應力剪切傳輸(SST)對迷宮節流件的流動特性進行了研究,五種湍流模型均采用穩態分析。計算方法選用Simple算法。

表4 邊界條件

2.3 網格敏感性分析

采用商業軟件Fluent進行數值模擬。在邊界層厚度參數相同的情況下,通過改變全局網格參數生成五種數量的網格。在鈉冷卻劑質量流量不變的情況下,分析了迷宮節流件壓降隨網格數的變化趨勢,如圖6所示,當總網格單元數達到786萬時,壓降增加幅度很小,計算網格的進一步細化并沒有改善數值結果。考慮到精度和收斂速度,選擇786萬個網格單元的計算網格作為工作網格。

圖6 網格敏感性分析

3 迷宮節流件的水力分析

3.1 節流件阻力系數

根據公式(4)可以求得節流件的阻力系數計算公式:

根據公式(8)和實驗數據可計算出阻力系數在兩種循環條件下隨雷諾數的變化。如圖7和圖8所示,在強迫循環條件下,阻力系數的變化很小,基本穩定在一個恒定的數值,即,。

圖8 自然循環下阻力系數隨雷諾數的變化

而在自然循環條件下,節流件的阻力系數會隨雷諾數的增加逐漸降低,并逐漸趨于平緩,對實驗數據進行回歸分析,可以得到阻力系數與雷諾數的關系式:

3.2 不同湍流模型的計算比較

3.2.1強制循環條件

圖9展示了使用不同湍流模型的計算結果與實驗結果之間的比較,在強制循環的高流量下,realizable-的計算結果與實驗結果最為吻合。標準-模型、S-A模型和SST模型的計算結果分別比實驗值高3.4%~5.42%、5.26%~7.56%和13.97%~18.45%。標準-模型是唯一計算結果低于實驗值的模型,比實驗值低6.37%~8.19%。因此,建議采用realizable-模型來模擬雷諾數244 323~352 796范圍內的迷宮式節流件。

同時當流量為 11.53 kg/s 的額定流量時,realizable-模型模擬的壓降為472.7 kPa,與實驗結果吻合良好,相對誤差僅為0.9%。驗證了高流量條件下數值模擬的可靠性。由于SFR的堆芯非常復雜,并非所有的結構設計都能通過實驗驗證。因此,數值模擬是對實驗研究的有效補充。

圖9 強迫循環條件下不同湍流模型的比較

影響湍流模型計算精度的除了網格數量,還有壁面法向網格距離,常常用無量綱數+表示。+定義為[15]:

——壁面距離。

在圖10中,采用realizable-湍流模型在+為10時的數值結果與實驗結果吻合良好,而+為30和5的計算結果的平均誤差則分別為-3.08%和-1.59%。結果表明,在強迫循環條件下,數值模擬對壁面法向網格距離非常敏感。

圖10 realizable k-ε模型采用不同y+值的計算結果比較

3.2.2自然循環條件

圖11為自然循環條件下不同湍流模型的計算結果,圖12是其各種湍流模型的計算誤差。如圖11和圖12所示,自然循環下五種湍流模型的表現與強迫循環下略有不同。首先,realizable-模型的計算精度降低;同時,標準-模型的精度有提高,誤差保持在5%以內;SST和標準-與實驗值仍有較大誤差。但隨著流量的減小,SST模型的計算精度有提高的趨勢,而標準-模型的計算精度進一步降低。這是因為隨著流速的降低,湍流的壁面粘性變得越來越突出,而在高流速下,壁面分子粘性并不明顯。因此流量的降低使得像SST等對近壁面剪切流動有更好模擬效果的湍流模型具有小的誤差。

圖11 自然循環條件下不同湍流模型的比較

圖12 自然循環條件下不同湍流模型的計算誤差

對于+的測試,采用低流量下表現最好的標準-湍流模型。如圖13所示,當流量大于0.4 kg/s時,+=30更合適,而當流量小于0.4 kg/s時,只有+=10的誤差可以在3%以內。

圖13 標準k-w模型采用不同y+值的計算結果比較

3.3 節流件流場特性分析

在迷宮型節流件中,液體鈉流依次被節流板阻擋,其中一部分可用壓頭被轉換為動能。隨后,動能被湍流渦旋消散并轉化為熱能[16]。為了直觀的分析流體速度場,創建了兩個通過軸線且相互垂直的截面,截面和截面。兩個剖面圖如圖14所示。

圖15顯示了額定流量為11.5kg/s下截面上的壓力和流速等值線。如圖15(a)所示,當鈉冷卻劑依次流過節流板時,壓力逐漸降低,流速周期性變化。每當液態鈉流經節流板頂部時,其壓力就會急劇下降,每個節流板頂部都會出現一定的低壓區域,但兩個節流板之間的間隙中壓力則變化不大。如圖15(b)所示,每個節流板間隙間的速度分布類似,均存在一個低流速區域,其中第一個間隙最為特殊,其低速區域明顯大于后續的節流板,同時第一個間隙間的最大速度也低于后續間隙中的最大速度。而在后續的間隙中,因為流道的180°轉向,液態鈉會沖擊下一個節流板的底部,最大流速便出現在沖擊處。在每個間隙的中間區域,流速會出現明顯分層現象。

圖14 XY截面和ZX截面

圖16顯示了從垂直于截面和截面的角度觀察到的三維流線。從圖16(a)中可以看出,在每個節流板間隙中,液態鈉將形成湍流漩渦。由于節流件流道復雜曲折,渦流不是簡單的縱向湍流渦流或橫向二次流,而是類似兩者的組合。如圖16(b)所示,每當液態鈉流經一個節流板的頂部時,由于流動截面的縮小和下一個節流板的阻擋,液態鈉向兩側擴散,形成橫向的漩渦。

圖15 XY界面上壓力和速度云圖

圖16 節流件流線圖

然而,結合圖16(a),在液態鈉橫向擴散時,它也也會形成縱向的漩渦。這兩股流動疊加在節流板間隙中,形成兩個傾斜的漩渦。這兩個漩渦便是形成間隙間低流速區域及流速分層的主要原因。

3.4 兩種循環條件下流場比較

圖17顯示了低流量0.29 kg/s和額定流量11.53 kg/s之間的流速等值線云圖比較。正如之前的討論,每個節流板頂部后面都會有回流區域。如圖17(a)所示,強迫循環條件下,第一個節流板后的回流較大,幾乎占據整個節流板的后側。從第二個節流板開始,回流區域迅速減少,并一直維持到最后一個節流板。而在自然循環條件下,盡管每個節流板后仍有回流,但這些區域不會很快收縮。如圖17(b)所示,在第二個和第三個節流板后面仍然有一個比較大的回流區域。同時從圖17(c)和圖17(d)還可以看出,低流量下的速度場更加無序,每個間隙之間的主流區域更小。這表明在低流速下,壁面對流動的擾動更為顯著。這也使得對于邊界層流動和低雷諾數流動具有良好模擬精度的標準-模型在自然循環下具有更好的精度。

圖17 流場速度云圖

4 結論

在本論文中,對迷宮節流件在自然和強制循環條件下進行了實驗和數值研究,得出以下結論。

(1)通過水力試驗,獲得節流件在自然循環條件下阻力系數與雷諾數的關系式,同時得到在強制循環條件下,節流件的阻力系數穩定在144.17左右。

(2)對節流件在兩種循環條件下進行數值模擬,通過與實驗結果的比較,推薦壁面法向網格間距+的取值為10。同時得出Realizable-湍流模型在強迫循環下具有最好的計算精度,而標準的-湍流模型在自然循環下具有最好的計算精度。

(3)節流板間隙內的傾斜湍流漩渦引起的能量耗散是節流件產生壓降的主要原因。同時,漩渦使每個節流板的在后側形成低速回流區域。回流區域會隨著流體不斷經過節流板而逐漸縮小。

(4)兩種循環條件下模擬結果的差異主要集中在二次回流區域大小的不同。在自然循環的低流量下,各節流板后的回流區域明顯大于高流量下的回流區域,且收縮速度較慢。

[1] USDOE.A Technology Roadmap for Generation IV Nuclear Energy Systems[R]. 2002:239-241.

[2] Fang X,Wang Y Z,Diao A N,et al.Structural Optimization Design of Labyrinth Seal Based on Software FLUENT[J]. Fluid Machinery,2013.

[3] Rhode D L,Sobolik S R.Simulation of Subsonic Flow Through a Generic Labyrinth Seal Cavity[C]. Asme International Gas Turbine Conference & Exhibit.American Society of Mechanical Engineers,1985.

[4] Jing Y,Luis S A.On the Influence of the Entrance Section on the Rotordynamic Performance of a Pump Seal with Uniform clearance:a Sharp Edge VS.a Round Inlet[J]. Journal of Engineering for Gas Turbines & Power,2018.

[5] Szymanski A,Wroblewski W,Fraczek D,et al. Optimization of the Straight-Through Labyrinth Seal With a Smooth Land[J]. Journal of Engineering for Gas Turbines & Power,2018,140(12):122503.1-122503.9.

[6] Zhang W,Wu K,Gu C,et al.Swirl Brakes Optimization for Rotordynamic Performance Improvement of Labyrinth Seals Using Computational Fluid Dynamics Method[J]. Tribology International,2021,159(2):106990.

[7] Hai-min,WANG,Yue,et al.Numerical simulation of flow characteristics for a labyrinth passage in a pressure valve[J]. 水動力學研究與進展:英文版,2016(4):8.

[8] Asok S P,Sankaranarayanasamy K,Sundararajan T,et al.Neural network and CFD-based optimisation of square cavity and curved cavity static labyrinth seals[J]. Tribology International,2007,40(7):1204-1216.

[9] Wang D,Shi C,Cai X.Numerical calculation of the influence of axial position on the performance of stepped labyrinth seals[C]. 6TH INTERNATIONAL CONFERENCE ON COMPUTER-AIDED DESIGN,MANUFACTURING,MODELING AND SIMULATION(CDMMS 2018),2018.

[10] Qin H,Lu D,Zhong D,et al.Experimental and numerical investigation for the geometrical parameters effect on the labyrinth-seal flow characteristics of fast reactor fuel assembly[J]. Annals of nuclear energy,2020,135(Jan.):106964.1-106964.14.

[11]宋怡,陸道綱,秦亥琦,等.迷宮型節流件壓降的數值模擬研究[J]. 核科學與工程,2020,40(2):9.

[12] Duenas Dobrowolski J,Gawliński,Marek,Paszkowski M,et al.Experimental Study of Lubricating Grease Flow inside the Gap of a Labyrinth Seal Using Microparticle Image Velocimetry[J]. Tribology Transactions,2016:1-10.

[13] Untaroiu A,Goyne C P,Untaroiu C D,et al.Computational Modeling and Experimental Investigation of Static Straight-Through Labyrinth Seals[C]. Asme International Mechanical Engineering Congress & Exposition,2008.

[14] Liu Y Z.Investigation on the flow and temperature fields of coolant in a fuel subassembly foe sodium cooled fast reactor[J]. China Institute of Atomic Energy,2007.

[15] Jing C,Zhang D,Ping S,et al.CFD investigation on thermal-hydraulic behaviors of a wire-wrapped fuel subassembly for sodium-cooled fast reactor[J]. Annals of Nuclear Energy,2018,113:256-269.

[16] Rhode D L,Sobolik S R.Simulation of Subsonic Flow Through a Generic Labyrinth Seal Cavity[C]. Asme International Gas Turbine Conference & Exhibit.American Society of Mechanical Engineers,1985.

Study on Hydraulic Characteristics of Labyrinth Throttle in Sodium Cooled Fast Reactor under Natural Circulation

ZHOU Zhiwei1,CAO Feng2,LU Daogang2,CAO Qiong2,*

(1. China institute of atomic energy,Beijing 102413,China;2. Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy,North China Electric Power University,Beijing 102206,China)

Labyrinth throttles are often installed in the flow bypass of the main vessel cooling system and pump cooling system of sodium cooled pool fast reactor, so that the flow entering these bypass reaches the design value, so as to effectively cool the equipment and ensure the safety of fast reactor. At present, there are few domestic research on the resistance characteristics of throttle under natural circulation.Therefore, this paper makes experimental and numerical research on labyrinth throttle under natural circulation and forced circulation. In the experiment, liquid water with constant temperature(84 ℃)was used as the working medium, and the resistance coefficients of throttling parts under two cycle conditions were obtained.Then, the throttle is numerically simulated, and the numerical simulation results under different turbulence models and the normal grid spacing of the wall are compared with the experimental results, which verifies the validity of the numerical simulation.By comparing the calculation results and flow field details of different turbulence models, a turbulence model suitable for the numerical simulation of labyrinth throttle under two conditions is obtained.At the same time, the difference of the flow field of the throttle under the two conditions is compared.It is found that the reflux vortex is the main energy consumption mode of the throttle, and the difference of the simulation results under the two conditions is mainly concentrated in the size of the secondary reflux area.Under natural circulation, the reflux area is significantly larger than that under forced circulation, and the contraction speed is slow.

Labyrinth throttle; Hydraulic characteristics; Turbulence model; Numerical simulation

TL48

A

0258-0918(2023)05-0979-10

2022-06-27

周志偉(1986—),男,湖南邵陽人,博士研究生,現從事核反應堆熱工水力學方面研究

曹 瓊,E-mail:Caoqiong@ncepu.edu.cn

猜你喜歡
實驗模型
一半模型
記一次有趣的實驗
微型實驗里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 最新亚洲人成无码网站欣赏网 | 福利视频一区| 精品国产网| 日韩亚洲综合在线| 国产精品偷伦在线观看| 国产在线视频自拍| 国产在线视频欧美亚综合| 亚瑟天堂久久一区二区影院| 日韩精品高清自在线| 91外围女在线观看| 国产精品黄色片| 91丝袜乱伦| 成年免费在线观看| 欧美亚洲国产精品久久蜜芽| 最新日本中文字幕| 人妻丰满熟妇av五码区| 2020亚洲精品无码| 一级在线毛片| 亚洲综合18p| 制服丝袜一区二区三区在线| 国产毛片一区| 亚洲三级片在线看| 漂亮人妻被中出中文字幕久久| 中文字幕在线观| 日韩麻豆小视频| 国产精品亚洲日韩AⅤ在线观看| 在线免费亚洲无码视频| 久爱午夜精品免费视频| 国产a v无码专区亚洲av| 99久久亚洲综合精品TS| 91人人妻人人做人人爽男同| 欧美另类精品一区二区三区| 成年人福利视频| 又黄又爽视频好爽视频| 99视频在线免费观看| 欧美精品不卡| 香蕉eeww99国产在线观看| 狠狠亚洲婷婷综合色香| 国产精品亚洲天堂| 精品无码日韩国产不卡av| 九九香蕉视频| 国产99视频在线| 国产精品深爱在线| 亚洲无码高清免费视频亚洲| 日韩精品成人在线| 亚洲中久无码永久在线观看软件| 手机精品福利在线观看| 九色在线观看视频| 91香蕉国产亚洲一二三区| 日本国产在线| 亚洲有无码中文网| 2019国产在线| 色网站免费在线观看| 免费无码又爽又黄又刺激网站| 午夜精品区| 思思热在线视频精品| 91国内在线视频| 午夜激情福利视频| 在线播放国产一区| 免费av一区二区三区在线| 久久福利网| 成人在线视频一区| 免费无码在线观看| 亚洲无码一区在线观看| 免费毛片网站在线观看| 欧美日韩免费观看| 国产免费观看av大片的网站| 欧美精品H在线播放| 呦女亚洲一区精品| 日韩精品一区二区三区swag| h视频在线观看网站| 亚洲永久精品ww47国产| 九色在线观看视频| 老熟妇喷水一区二区三区| 国产精品对白刺激| 天天做天天爱夜夜爽毛片毛片| 91偷拍一区| 第一区免费在线观看| 一级爱做片免费观看久久| 日韩亚洲综合在线| 精品国产美女福到在线不卡f| 欧类av怡春院|