亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        鈉冷快堆繞絲環(huán)形燃料組件流動特性分析

        2023-12-16 05:38:28袁顯寶雷蔭先郭躍峰周建軍杜曉超張永紅張彬航楊森權(quán)
        核科學(xué)與工程 2023年5期
        關(guān)鍵詞:模型研究

        袁顯寶,雷蔭先,郭躍峰,*,周建軍,杜曉超,張永紅,張彬航,譚 超,楊森權(quán)

        鈉冷快堆繞絲環(huán)形燃料組件流動特性分析

        袁顯寶1,雷蔭先2,3,郭躍峰2,3,*,周建軍2,3,杜曉超2,3,張永紅2,3,張彬航2,3,譚超4,楊森權(quán)4

        (1. 三峽大學(xué) 理學(xué)院,湖北 宜昌 443002;2. 三峽大學(xué) 機械與動力學(xué)院,湖北 宜昌 443002;3. 湖北省水電機械設(shè)備設(shè)計與維護重點實驗室 三峽大學(xué),湖北 宜昌 443002;4. 中核武漢核電運行技術(shù)股份有限公司,湖北 武漢 430074)

        鈉冷快堆通過金屬繞絲定位,繞絲在加強冷卻劑交混、促進對流換熱、減小機械振動等方面具有重要的作用。本文通過計算流體力學(xué)(CFD)建立鈉冷快堆環(huán)形燃料組件的相關(guān)數(shù)學(xué)物理模型,就繞絲所引起冷卻劑在子通道的流動特性展開研究。研究表明:冷卻劑在子通道內(nèi)的流動呈現(xiàn)周期性,周期為一個螺距(100 mm);中心子通道內(nèi)冷卻劑的流動受相鄰棒束的三個繞絲的影響,繞絲的交混作用更強,流動更加穩(wěn)定;子通道橫向流最大為1.74 m/s,橫向流最大的位置隨著繞絲螺旋結(jié)構(gòu)纏繞位置的不同而變化;子通道間隙內(nèi)的歸一化橫向流具有周期性。本文為鈉冷快堆環(huán)形燃料組件的進一步優(yōu)化設(shè)計提供數(shù)據(jù)支撐,具有一定的研究意義。

        鈉冷快堆;繞絲;橫向流;計算流體力學(xué)

        鈉冷快堆作為第四代堆型之一,針對其堆芯熱工水力方面的研究具有重要的意義。鈉冷快堆燃料棒采用金屬繞絲定位,繞絲在充分交混冷卻劑、加強對流換熱等方面具有重要的作用。因此,關(guān)于繞絲的研究受到眾多學(xué)者的青睞。

        Novendstern[1]、Rehme[2]、Engel等[3]、Cheng and Todreas等[4]通過實驗研究不同棒束組件的冷卻劑壓降的變化,并根據(jù)實驗擬合得到雷諾數(shù)與摩擦因子的關(guān)系式。Liang等[5]通過實驗研究37棒束帶繞絲燃料組件壓降及流動分布特性,并與經(jīng)驗關(guān)系式進行對比。叢騰龍等[6]對鉛鉍快堆不同繞絲數(shù)量時燃料組件橫向流特性展開研究。姜文殊等[7]基于CFD數(shù)值模擬,研究不同子通道的平均速度壓降等。何明翰等[8]以CiADS為研究對象,通過對比不同網(wǎng)格下棒束通道的橫向流特性,改進橫流特性經(jīng)驗公式。王婧婕[9]、范大軍[10]等均對子通道橫向流展開研究。在上述關(guān)于繞絲流動特性的實驗與數(shù)值模擬主要集中在棒狀燃料組件方面,未涉及環(huán)形燃料組件。環(huán)形燃料組件具有內(nèi)外兩個冷卻劑流道,可充分帶走燃料芯塊產(chǎn)生的熱量,降低冷卻劑溫度,提升堆芯功率密度。美國[11,12]、韓國[13]、中國[14]等均對壓水堆使用環(huán)形燃料組件進行研究,證明其功率提升20%~50%是可行的??紤]到環(huán)形燃料的諸多優(yōu)點,試想可否將環(huán)形燃料應(yīng)用于第四代堆型。季松濤等[15]提出將環(huán)形燃料應(yīng)用于鈉冷快堆。為此,本文在鈉冷快堆環(huán)形燃料組件的基礎(chǔ)上,針對繞絲燃料組件的流動特性展開研究。

        本文通過計算流體力學(xué)(CFD)建立鈉冷快堆環(huán)形燃料組件的相關(guān)數(shù)學(xué)物理模型,就繞絲所引起的冷卻劑在子通道中的流動特性展開研究,涉及軸向流動及橫向流特性,為鈉冷快堆環(huán)形燃料組件的進一步研究,提供了參考,具有一定的研究意義。

        1 數(shù)值模型及設(shè)置

        7 棒束鈉冷快堆環(huán)形燃料組件的數(shù)值模型主要包括幾何模型、網(wǎng)格劃分、控制方程及湍流模型的選取、邊界條件的設(shè)置4個部分。

        1.1 幾何模型

        本文依據(jù)燃料組件尺寸不變、冷卻劑——核燃料體積比不變的原則設(shè)計了中國實驗快堆環(huán)形燃料組件。環(huán)形燃料組件的幾何參數(shù)如表1所示。選取中心的七棒束作為研究對象,繞絲均順時針纏繞在環(huán)形燃料棒上,繞絲出口纏繞角設(shè)置為30°,幾何模型截面圖如圖1所示。

        表1 環(huán)形燃料組件參數(shù)

        圖1 幾何模型截面示意圖

        1.2 網(wǎng)格劃分及無關(guān)性分析

        本文利用三維設(shè)計軟件Solidworks對繞絲環(huán)形燃料組件進行建模,使用ANSYS DM對幾何部分進行處理??紤]到繞絲與環(huán)形燃料棒之間的接觸為線接觸,二者接觸的位置易出現(xiàn)尖角,導(dǎo)致網(wǎng)格質(zhì)量降低。因此,為了簡化計算,將繞絲與環(huán)形燃料棒相向移動0.05 mm,使線接觸變?yōu)槊娼佑|,這樣微小的改變對計算結(jié)果的影響忽略不計[16]。采用Fluent Meshing繪制多面體網(wǎng)格,對于繞絲和子通道,由于其幾何復(fù)雜,形狀不規(guī)則,通過局部網(wǎng)格加密的方式生成多面體網(wǎng)格。在冷卻劑與環(huán)形燃料組件的接觸面處、冷卻劑與燃料盒壁面的接觸面處均設(shè)置了邊界層,按照增強型壁面函數(shù)的值來設(shè)定的第一層網(wǎng)格高度。其局部網(wǎng)格如圖2所示。

        圖2 局部網(wǎng)格示意圖

        針對本文幾何模型,選取不同幾何尺寸的四個例子,得到不同數(shù)量的網(wǎng)格,對比隨著軸向高度增加時冷卻劑溫度的變化,如圖3所示。隨著/m的增加,冷卻劑溫度逐漸增加。四個例子對應(yīng)的冷卻劑最高溫度的最大誤差僅為0.005%,所有的誤差均在可接受的范圍內(nèi)。從計算精度考慮,四個例子的精度均滿足計算要求;從計算時間考慮,選擇網(wǎng)格數(shù)量較少的例子節(jié)約計算時間。最終,選取例子1作為本文數(shù)值計算的網(wǎng)格模型,其網(wǎng)格數(shù)量為231萬。

        圖3 網(wǎng)格無關(guān)性分析

        1.3 控制方程及湍流模型

        通過求解控制方程得到七棒束繞絲環(huán)形燃料組件的流動特性,三大守恒方程如下:

        連續(xù)性方程:

        式中:——流體的密度;

        ,,——流體在,,三個方向上的速度。

        動量方程:

        式中:——流體黏性應(yīng)力;

        ——流體動力黏度;

        ——壓力。

        能量方程:

        式中:——流體溫度;

        ——流體導(dǎo)熱系數(shù);

        C——流體比熱容;

        t——源項。

        本文選取標(biāo)準(zhǔn)-模型作為湍流模型,該模型由Launder和Spalding[17]提出,包括湍動能和耗散率兩個輸運方程,其方程如下:

        (5)

        式中:——湍動能;

        ——湍流耗散系數(shù);

        t——流體的湍流黏度;

        、S——流動的變形率張量;

        其余參數(shù)均為常數(shù),分別為:C1.44,C1.92,μ0.09,k1.0,ε1.3。t、S關(guān)系式如下。

        1.4 邊界條件

        中國實驗快堆[18]采用液態(tài)金屬鈉作為冷卻劑,本文將鈉冷快堆冷卻劑物性設(shè)置為變物性,考慮鈉的密度、導(dǎo)熱系數(shù)、比定壓熱容p、黏度,其關(guān)系式如表2所示,以UDF的形式導(dǎo)入Fluent。

        普朗特數(shù)作為表征流體傳熱特性重要的無量綱參數(shù),較水和空氣,液態(tài)金屬鈉普朗特數(shù)遠(yuǎn)小于1,表現(xiàn)為導(dǎo)熱能力遠(yuǎn)大于對流擴散能力,溫度邊界層厚度遠(yuǎn)大于流動邊界層厚度,因此,F(xiàn)luent默認(rèn)的t0.85已不適用于液態(tài)金屬鈉。本文基于Subbotin[19]實驗,驗證Jischa[20]提出的湍流普朗特數(shù)t模型的適用性,如圖4所示。隨著貝克萊數(shù)的增加,努塞爾數(shù)逐漸增加。其中t0.85模型與Subbotin實驗誤差較大;Jischa模型與Subbotin實驗誤差較小。為此,本文選取Jischa模型作為t模型,并以UDF的形式導(dǎo)入。

        表2 冷卻劑物性參數(shù)

        根據(jù)中國實驗快堆的設(shè)計準(zhǔn)則[18],冷卻劑入口溫度為633.15 K,入口采用質(zhì)量流量——溫度入口,內(nèi)流場冷卻劑入口質(zhì)量流量為0.07 kg/s,子通道冷卻劑入口質(zhì)量流量為1 kg/s;出口設(shè)置為壓力出口,出口靜壓為0 Pa;環(huán)形燃料棒熱源設(shè)置為平均表面熱流密度,其值為1.49×106W/m2;繞絲、燃料盒壁面均為無滑移絕熱壁面。

        圖4 Prt模型驗證

        2 模型驗證

        流動特性作為研究燃料組件熱工水力特性的重要部分,從20世紀(jì)70年代開始,各國學(xué)者通過實驗研究子通道壓降特性得到眾多摩擦因子的經(jīng)驗關(guān)系式,如Rehme、Cheng and Todreas等。以上經(jīng)驗關(guān)系式主要針對棒狀燃料組件,為此,僅驗證本模型子通道部分,不考慮內(nèi)流場。根據(jù)幾何模型的/,/及雷諾數(shù)的范圍,選取與之相符合的關(guān)系式進行對比驗證,如圖5所示。隨著雷諾數(shù)的增加,摩擦因子逐漸減小。經(jīng)驗證:本文模型與修正后的Rehme關(guān)系式吻合較好,數(shù)值計算的結(jié)果具有一定的說服力。

        圖5 摩擦因子關(guān)系式驗證

        3 流動特性分析

        3.1 軸向流動特性

        圖6所示為子通道編號,1~6號子通道內(nèi)冷卻劑的流動主要受周圍3根燃料棒的影響,稱為中心子通道;7~12號子通道位于燃料盒壁面邊的位置,稱為邊子通道;13~18號子通道位于燃料盒角的位置,稱為角子通道。本文選取最有代表性的三個子通道,編號分別為:4、7、18。

        圖6 子通道編號

        各子通道冷卻劑沿軸向高度的平均速度變化如圖7所示。從整體上看:隨著軸向高度的增加,子通道冷卻劑的平均速度呈現(xiàn)周期性的變化。各子通道冷卻劑平均速度的變化曲線均出現(xiàn)較明顯的五個“峰”且周期為一個螺距,表現(xiàn)在圖7為/2,其值為12.5 mm。中心子通道內(nèi)的冷卻劑流動由于受到相鄰的三個燃料棒繞絲的影響,平均速度波動更加明顯,在一個周期內(nèi)出現(xiàn)三次小幅波動。

        從局部來看:位于邊子通道的冷卻劑平均速度變化幅度為2.37 m/s,遠(yuǎn)大于位于角子通道和中心子通道的冷卻劑平均速度的變化幅度,中心子通道的冷卻劑平均速度變化幅度最小,僅為1.12 m/s。就最大平均速度來看:邊子通道最大,為6.62 m/s;角子通道次之,為6.15 m/s;中心子通道內(nèi)最小,僅為6.09 m/s。由此可知:中心子通道內(nèi)的冷卻劑流動較其他子通道更加穩(wěn)定,冷卻劑流速波動較小。中心子通道內(nèi)冷卻劑的流動受相鄰棒束的三個繞絲的影響,繞絲的交混作用更強,使流動更加穩(wěn)定。

        圖7 各子通道冷卻劑平均速度

        3.2 橫向流特性分析

        圖8所示為112.5 mm、225 mm、337.5 mm、450 mm時冷卻劑橫向流云圖。由圖8可知,子通道橫向流最大為1.74 m/s,最小為0.17 m/s,相差一個數(shù)量級。橫向流最大的位置隨著繞絲螺旋結(jié)構(gòu)纏繞位置的不同而變化,表現(xiàn)為不同截面橫向流最大的位置具有差異性。對比不同截面的橫向流云圖,不難發(fā)現(xiàn):角子通道與燃料棒之間的狹縫處以及燃料棒與繞絲相接的位置橫向流最大。當(dāng)冷卻劑流經(jīng)角子通道與燃料棒之間的狹縫時,流通面積減小,流速增加,繞絲的交混作用使速度在軸、軸的分量增加,所以橫向流增加。當(dāng)冷卻劑流經(jīng)燃料棒繞絲位置時,繞絲的螺旋結(jié)構(gòu)使冷卻劑流向發(fā)生變化,繞絲位置處冷卻劑的流動出現(xiàn)“漩渦”,如圖9橫向流矢量圖所示。這種“漩渦”帶動冷卻劑流速增加,橫向流增加。

        圖8 橫向流

        為了更加深入地研究橫向流特性,本文選取不同位置的子通道間隙1#、2#、3#,如圖10所示。子通道間隙1#表示內(nèi)層燃料棒與外層燃料棒間隙;子通道間隙2#表示同一層燃料棒之間的間隙;子通道3#表示燃料棒與燃料盒壁面之間的間隙。在提取間隙1#、2#、3#的橫向流數(shù)據(jù)時,定義了歸一化橫向流的公式如式(8)所示。

        式中:——子通道間隙1#、2#、3#的寬度,該值會隨繞絲位置的變化而變化;

        圖11所示為一個螺距(100 mm)內(nèi),子通道間隙1#、2#、3#的歸一化橫向流變化曲線。在一個螺距內(nèi),繞絲會兩次繞過1#、2#、3#使間隙的寬度減小,歸一化橫向流曲線會有兩個位置接近于0,且兩個位置間隔半個螺距(50 mm),即歸一化橫向流變化呈現(xiàn)周期性。繞絲所占據(jù)角度為12°,所以歸一化橫向流為0的位置占據(jù)的角度為12°。在間隙1#時,繞絲分別在150°、330°繞過該間隙,對應(yīng)1#曲線中歸一化橫向流為0的位置;間隙2#時,繞絲分別在30°、210°繞過該間隙,對應(yīng)2#曲線中歸一化橫向流為0的位置;間隙3#時,繞絲分別在60°、240°繞過該間隙,對應(yīng)3#曲線中歸一化橫向流為0的位置。

        圖9 橫向流矢量圖

        圖10 子通道間隙命名示意圖

        圖11 歸一化橫向流

        4 結(jié)論

        本文通過計算流體力學(xué)(CFD)建立鈉冷快堆環(huán)形燃料組件的相關(guān)數(shù)學(xué)物理模型,就繞絲所引起的冷卻劑在子通道中的流動特性展開研究,結(jié)論如下:

        (1)冷卻劑在子通道內(nèi)的流動呈現(xiàn)周期性,周期為一個螺距(100 mm)。中心子通道內(nèi)冷卻劑的流動受相鄰棒束三個繞絲的影響,繞絲的交混作用更強,流動更加穩(wěn)定。

        (2)子通道橫向流最大為1.74 m/s,且橫向流最大的位置隨著繞絲螺旋結(jié)構(gòu)纏繞位置的不同而變化。當(dāng)冷卻劑流經(jīng)燃料棒繞絲位置時,繞絲位置處冷卻劑的流動出現(xiàn)“漩渦”使冷卻劑流速增加,橫向流增加。

        (3)子通道間隙內(nèi)的歸一化橫向流具有周期性。本文為鈉冷快堆環(huán)形燃料組件的下一步優(yōu)化,提供了參考,具有一定的研究意義。

        [1] Novendstern E H.Turbulent flow pressure drop model for fuel rod assemblies utilizing a helical wire-wrap spacer system[J]. Nuclear Engineering and Design,1972,22(1):28-42.

        [2] Rehme,Klaus.Pressure Drop Correlations for Fuel Element Spacers[J]. Nuclear Technology,1973,17(1):15-23.

        [3] Engel F C,Markley R A,Bishop A A.Laminar,Transition,and Turbulent Parallel Flow Pressure Drop Across Wire-Wrap-Spaced Rod Bundles[J]. Nuclear science and engineering:the journal of the American Nuclear Society,1979,69(2):290-296.

        [4] Cheng S K,Todreas N E.Hydrodynamic models and correlations for bare and wire-wrapped hexagonal rod bundles—bundle friction factors,subchannel friction factors and mixing parameters[J]. Nuclear engineering and design,1986,92(2):227-251.

        [5] Liang Y,Zhang D,Chen Y,et al.An experiment study of pressure drop and flow distribution in subchannels of a 37-pin wire-wrapped rod bundle[J]. Applied Thermal Engineering,174.

        [6] 叢騰龍,王俊杰,肖瑤,等. 鉛鉍冷卻繞絲燃料組件橫流特性分析[J]. 原子能科學(xué)技術(shù),2022,56(12):2725-2734.

        [7] 姜文殊,梁峻銘,張會勇,等. 液態(tài)鉛鉍合金在繞絲燃料棒組件中的水力特性數(shù)值模擬[J]. 東北電力大學(xué)學(xué)報,2022,42(03):74-82101.

        [8] 何明翰,范大軍,李榮杰,等. CiADS帶繞絲燃料棒束通道橫流特性[J]. 原子核物理評論,2022,39(02):258-265.

        [9] 王婧婕.液態(tài)鉛鉍合金在繞絲燃料棒組件中的流動與傳熱機理研究[D]. 北京:北京化工大學(xué),2021.

        [10]范大軍.鉛冷快堆帶繞絲燃料棒束通道流動特性研究[D]. 蘭州中國科學(xué)院大學(xué)(中國科學(xué)院近代物理研究所),2021.

        [11] Lahoda E,Mazzoccoli J,Beccherle J.High-power-density annular fuel for pressurized water reactors:manufacturing costs and economic benefits[J]. Nuclear technology,2007,160(1):112-134.

        [12] Feng D,Morra P,Sundaram R,et al.Safety analysis of high-power-density annular fuel for PWRs[J]. Nuclear technology,2007,160(1):45-62.

        [13] Shin C H,Chun T H,Oh D S,et al.Thermal hydraulic performance assessment of dual-cooled annular nuclear fuel for OPR-1000[J]. Nuclear Engineering and Design,2012,243:291-300.

        [14]馮海寧,趙瑞瑞,王虹,等. 壓水堆環(huán)形燃料組件研發(fā)綜述[J]. 中國核電,2020,13(06):759-764.

        [15]季松濤,韓智杰,何曉軍,等. 壓水堆環(huán)形燃料組件研發(fā)進展[J]. 原子能科學(xué)技術(shù),2020,54(S1):240-245.

        [16] Merzari E,Pointer W D,Smith J G,et al. Numerical simulation of the flow in wire-wrapped pin bundles:Effect of pin-wire contact modeling[J]. Nuclear engineering and design,2012,253:374-386.

        [17] Launder B E,Spalding D B.The numerical computation of turbulent flows[M]. Numerical prediction of flow,heat transfer,turbulence and combustion.Pergamon,1983:96-116.

        [18]徐銤.快堆熱工流體力學(xué)[M]. 北京:原子能出版社,2011.

        [19] Subbotin V I,Papovyants A K,Kirillov P L,et al.A study of heat transfer to molten sodium in tubes[J]. Soviet Atomic Energy,1963,13(4):991-994.

        [20] Jischa M,Rieke H B.About the prediction of turbulent Prandtl and Schmidt numbers from modeled transport equations[J]. International Journal of Heat and Mass Transfer,1979,22(11):1547-1555.

        Flow Characteristic Analysis of the Wire-wrapped Annular Fuel Assembly for the Sodium-Cooled Fast Reactor

        YUAN1Xianbao,LEI Yinxian2, 3,GUO Yuefeng2, 3, *,ZHOU Jianjun2, 3,DU Xiaochao2, 3,ZHANG Yonghong2, 3,ZHANG Binhang2, 3,TAN Chao4,YANG Senquan4

        (1.School of Science,China Three Gorges University,Yichang of Hubei Prov. 443002,China;2. School of Mechanical Engineering,China Three Gorges University,Yichang of Hubei Prov. 443002,China;3. Hubei Key Laboratory of Hydroelectric Machinery Design & Maintenance,Yichang of Hubei Prov. 443002,China 4. China Nuclear Power Operation Technology Corporation,Ltd.,Wuhan of Hubei Prov. 430074,China)

        The sodium-cooled fast reactor is positioned by metal wires, which plays an important role in enhancing coolant mixing, promoting convection heat transfer and reducing mechanical vibration. In this study, the mathematical and physical model of the annular fuel assembly of the sodium-cooled fast reactor is established by computational fluid dynamics (CFD), and the coolant flow characteristics in the sub-channels caused by wire-wrap are studied. The results show that the coolant flow in the sub-channels is periodic, and the period is one pitch (100 mm). The coolant flow in the central sub-channel is affected by the three wires of adjacent rod bundles, and the mixing effect of the wires is stronger and the flow is more stable. The maximum transverse flow is 1.74 m/s, and the position of the maximum transverse flow varies with the winding position of the spiral structure. This paper provides a reference for further optimization design of annular fuel assembly for the sodium-cooled fast reactor, and has certain research significance.

        Sodium-cooled fast reactor; Wire-wrap; Transverse flow; Computational fluid dynamics

        TL48

        A

        0258-0918(2023)05-1150-08

        2022-10-14

        國家自然科學(xué)基金項目(12175116)

        袁顯寶(1974—),男,湖北興山人,教授,博士,現(xiàn)從事反應(yīng)堆熱工物理方面研究

        郭躍峰,E-mail:verona1206@163.com

        猜你喜歡
        模型研究
        一半模型
        FMS與YBT相關(guān)性的實證研究
        2020年國內(nèi)翻譯研究述評
        遼代千人邑研究述論
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        視錯覺在平面設(shè)計中的應(yīng)用與研究
        科技傳播(2019年22期)2020-01-14 03:06:54
        EMA伺服控制系統(tǒng)研究
        新版C-NCAP側(cè)面碰撞假人損傷研究
        3D打印中的模型分割與打包
        熟女俱乐部五十路二区av| 阴唇两边有点白是怎么回事| 嫩草伊人久久精品少妇av| 人妻聚色窝窝人体www一区| 久久久久99精品成人片试看| 国产高清一级毛片在线看| 极品少妇高潮在线观看| 亚洲色精品三区二区一区| 射死你天天日| 亚洲ⅤA中文字幕无码| 亚洲av毛片在线播放| 日韩人妻熟女中文字幕a美景之屋 国产suv精品一区二区四 | 免费人妻精品一区二区三区| 狠狠色丁香婷婷综合潮喷| 无码人妻丰满熟妇啪啪网站| 老少交欧美另类| 久久久调教亚洲| 蜜桃免费一区二区三区| 伊人久久大香线蕉av不卡| 天堂…在线最新版资源| 亚洲AⅤ无码国精品中文字慕| 99国语激情对白在线观看| 色婷婷一区二区三区久久亚洲| 久久夜色国产精品噜噜亚洲av| 日本熟妇hdsex视频| 国产男女猛烈无遮挡免费视频| 国产白浆精品一区二区三区| 亚洲av午夜一区二区三| 精品国产一区二区三区免费| 中文字幕一区二区三区四区在线| 伊人久久亚洲综合av影院| 18国产精品白浆在线观看免费| 野外性史欧美k8播放| 日韩不卡无码三区| 青青河边草免费在线看的视频| 亚洲国产天堂一区二区三区| 97人妻视频妓女网| 成人影院视频在线播放| 国产日韩精品欧美一区喷水| 亚洲成aⅴ人在线观看 | 日韩中文网|