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

快堆帶繞絲棒束組件低雷諾數(shù)下的水力特性分析

2019-08-29 03:04:24程道喜齊曉光杜開文翟偉明
原子能科學(xué)技術(shù) 2019年8期
關(guān)鍵詞:模型

程道喜,齊曉光,杜開文,翟偉明

(中國(guó)原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部,北京 102413)

國(guó)際上的快堆組件大多采用帶繞絲棒束的形式,燃料棒可布置得很緊湊。快堆燃料組件棒束采用帶繞絲的結(jié)構(gòu)定位的同時(shí),可加強(qiáng)各子通道間冷卻劑的交混,強(qiáng)化燃料棒與冷卻劑的換熱,但繞絲的存在會(huì)增大棒束段的阻力。各國(guó)如美國(guó)、法國(guó)、日本、意大利等均針對(duì)燃料組件的水力特性(主要是棒束段的水力特性)進(jìn)行了大量的實(shí)驗(yàn)研究,并提出了一系列帶繞絲棒束組件的阻力壓降計(jì)算模型。Engel等[1]、Cheng和Todreas[2]利用實(shí)驗(yàn)總結(jié)了壓降經(jīng)驗(yàn)關(guān)系式。此外劉一哲等[3]在CRT模型和Engel等研究的基礎(chǔ)上,針對(duì)中國(guó)實(shí)驗(yàn)快堆(CEFR)燃料組件的具體情況提出了ICRT壓降關(guān)系式,用以計(jì)算冷卻劑在湍流區(qū)、過渡流區(qū)和層流區(qū)的棒束壓降。

近年來,日益強(qiáng)大的計(jì)算能力使計(jì)算流體力學(xué)(CFD)手段得以應(yīng)用到堆芯組件流動(dòng)及傳熱分析的相關(guān)領(lǐng)域,如壓水堆組件格架熱工水力特性研究[4]。由于快堆燃料組件結(jié)構(gòu)復(fù)雜,對(duì)于快堆燃料組件內(nèi)部流場(chǎng)及溫度場(chǎng)的CFD研究起步較晚。Ahmad等[5-14]利用CFD方法從不同角度進(jìn)行了快堆帶繞絲棒束組件的研究。綜合各項(xiàng)研究結(jié)果,近些年國(guó)外學(xué)者利用CFD程序?qū)於呀M件帶繞絲棒束的三維流動(dòng)特性以及傳熱特性進(jìn)行了大量的研究,但研究大多集中在湍流區(qū),而對(duì)于處于過渡區(qū)甚至是層流區(qū)的研究很少。

層流及過渡流動(dòng)這種低雷諾數(shù)(Re)流動(dòng)工況對(duì)于快堆的安全運(yùn)行是非常重要的。不僅在正常停堆或事故停堆后組件內(nèi)部會(huì)出現(xiàn)層流及過渡流動(dòng)的情況,甚至在反應(yīng)堆正常運(yùn)行時(shí)的某些組件內(nèi)也會(huì)出現(xiàn)。本文將利用CFX計(jì)算程序?qū)EFR帶繞絲棒束組件低Re下的水力特性進(jìn)行分析。

1 幾何模型

快堆帶繞絲棒束組件內(nèi)的燃料棒通常按照三角形柵格的形式排列,燃料棒之間利用金屬繞絲定位。繞絲按照一定的螺距螺旋式地順時(shí)針纏繞在燃料棒上。繞絲的存在大幅增加了幾何建模的難度。在實(shí)際結(jié)構(gòu)中,繞絲和燃料棒的接觸為相切接觸,在幾何建模時(shí)需將這個(gè)點(diǎn)接觸變成面接觸(即將繞絲與燃料棒的接觸由相切變?yōu)橄嘟?,以便進(jìn)行計(jì)算流體域模型建立及網(wǎng)格劃分。

CEFR燃料組件由61根帶繞絲的燃料棒組成,其中組件棒束段的長(zhǎng)度為1 350 mm,對(duì)整盒組件進(jìn)行完全模擬需要巨大的計(jì)算資源。在有限條件的情況下,本文選用較少螺距的7、19以及61根帶繞絲的燃料棒組成的組件作為對(duì)象來進(jìn)行模擬計(jì)算。圖1為1個(gè)螺距長(zhǎng)度的61根帶繞絲燃料棒的組件,其燃料棒和繞絲的幾何尺寸與CEFR燃料組件的一致:燃料棒直徑Dr=6 mm,繞絲直徑Dw=0.95 mm,燃料棒中心距Pr=7 mm,螺距H=100 mm。

圖1 計(jì)算模型示意圖Fig.1 Diagram of calculation model

2 網(wǎng)格劃分

考慮到快堆帶繞絲棒束組件結(jié)構(gòu)的復(fù)雜性,不可能采用結(jié)構(gòu)化網(wǎng)格進(jìn)行網(wǎng)格劃分,本文采用幾何適應(yīng)性好的非結(jié)構(gòu)化網(wǎng)格,并對(duì)繞絲周圍進(jìn)行加密。

繞絲的存在使組件結(jié)構(gòu)尺寸的尺度變化大,繞絲與相鄰燃料棒之間的間隙非常小,這樣造成網(wǎng)格劃分時(shí)非結(jié)構(gòu)化網(wǎng)格的數(shù)量非常龐大,在有限計(jì)算資源的情況下很難對(duì)組件整體進(jìn)行計(jì)算分析,僅能對(duì)很少螺距的帶繞絲棒束組件進(jìn)行模擬計(jì)算。為能獲得組件充分發(fā)展段的水力特性,在僅能對(duì)很少螺距的組件進(jìn)行模擬計(jì)算的情況下如果采用常規(guī)進(jìn)出口條件進(jìn)行計(jì)算,很可能由于入口段的影響,組件內(nèi)的流動(dòng)達(dá)不到充分發(fā)展。因此需選用平移周期性邊界條件來進(jìn)行處理,利用前一次迭代的出口的結(jié)果作為下一次迭代的入口條件,直至達(dá)到穩(wěn)定收斂時(shí)進(jìn)出口的結(jié)果一致以模擬組件內(nèi)充分發(fā)展的流動(dòng)。

3 邊界條件

3.1 1個(gè)螺距帶繞絲棒束組件

為得到充分發(fā)展段的帶繞絲棒束組件的水力特性,結(jié)合繞絲的周期性,在計(jì)算中對(duì)進(jìn)出口平面采用平移周期邊界條件,指定通過的質(zhì)量流量。由于采用了周期邊界條件,僅需計(jì)算1個(gè)螺距的帶繞絲棒束組件即可,可在有限的計(jì)算資源下計(jì)算更多數(shù)量的燃料棒組成的組件。

計(jì)算中使用液態(tài)金屬鈉作為工質(zhì),CEFR燃料組件入口鈉溫度為360 ℃,出口鈉溫度為530 ℃,故選取組件出入口平均溫度445 ℃為定性溫度確定鈉的物性。

3.2 4個(gè)螺距帶繞絲棒束組件

通過對(duì)4個(gè)螺距帶繞絲棒束組件的水力特性進(jìn)行分析,可得到其入口、出口段的長(zhǎng)度以及充分發(fā)展段水力特性變化規(guī)律。4個(gè)螺距的帶繞絲棒束組件需要大量的計(jì)算資源,本文僅進(jìn)行了7根燃料棒組成的組件的水力特性分析。

由于在此部分計(jì)算中要對(duì)組件出入口情況進(jìn)行分析,故不采用周期性邊界條件。設(shè)置入口質(zhì)量流量,同時(shí)出口為開口邊界,邊界壓力設(shè)為0,其他條件與3.1節(jié)中一致。

4 水力特性結(jié)果

4.1 1個(gè)螺距帶繞絲棒束組件

由于繞絲的交混作用,帶繞絲棒束組件在很低Re時(shí)已進(jìn)入過渡工況。本文分別使用層流模型和SST轉(zhuǎn)捩模型進(jìn)行計(jì)算并將計(jì)算結(jié)果與國(guó)際上公開發(fā)表的經(jīng)驗(yàn)Engel關(guān)系式[1]、Cheng關(guān)系式[2]以及劉一哲針對(duì)CEFR組件提出的ICRT關(guān)系式[3]的計(jì)算結(jié)果進(jìn)行比較,從而得到與實(shí)際情況比較接近的計(jì)算模型。

Engel關(guān)系式如下。

層流:

(1)

湍流:

(2)

過渡流動(dòng):

400≤Re≤5 000

(3)

式中:f為摩擦阻力系數(shù);Ψ為間斷因子,是一連接湍流區(qū)和層流區(qū)的參數(shù),Ψ=(Re-400)/4 600。

Cheng關(guān)系式如下。

層流:

(4)

湍流:

(5)

過渡流動(dòng):

ReL

(6)

其中:

式中:CfL、CfT分別為層流和湍流關(guān)系式的系數(shù);ReL為層流臨界雷諾數(shù);ReT為湍流臨界雷諾數(shù);Pt為棒束柵格距。

ICRT關(guān)系式如下。

層流:

(7)

湍流:

f=g1d1fs1M1X1+g2d2fs2M2X2

Re>5 000

(8)

過渡流動(dòng):

400≤Re≤5 000

(9)

式中:φ為與水力當(dāng)量粗糙度相關(guān)的參數(shù);g1、g2為面積分配因子;fs1、fs2為無繞絲通道的表面摩擦因子;d1、d2為棒束通道平均水力直徑與子通道水力直徑的比值;M1、M2為由繞絲引起的摩擦倍數(shù);X1、X2為不同通道的割流參數(shù)[3]。對(duì)于CEFR,φ=25.8,M1=1.124,湍流區(qū)參數(shù)本文不做詳細(xì)介紹。

考慮到各類型組件存在不同的棒束數(shù)量,本文計(jì)算了61、19及7根燃料棒的組件。

1) 61根帶繞絲棒束組件的水力特性

圖2 61根帶繞絲棒束組件水力特性計(jì)算結(jié)果與實(shí)驗(yàn)以及各關(guān)系式的比較Fig.2 Comparison of hydraulic characteristic calculation result of 61 wire-wrapped rod bundle assembly with experimental result and different relational expressions

61根帶繞絲棒束組件水力特性計(jì)算結(jié)果如圖2所示。從圖2可看出,利用SST轉(zhuǎn)捩模型計(jì)算的結(jié)果在大多情況下與俄羅斯實(shí)驗(yàn)結(jié)果[15]符合較好,尤其是在Re為300~2 000時(shí),計(jì)算值與實(shí)驗(yàn)值的偏差較小,隨著Re變小或Re變大時(shí)這一偏差逐步增大。層流模型計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)以及經(jīng)驗(yàn)關(guān)系式的偏差均較大;而在上述各經(jīng)驗(yàn)關(guān)系式中,Engel關(guān)系式的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果符合最好,在絕大多數(shù)情況下與實(shí)驗(yàn)值的偏差小于5%;而ICRT關(guān)系式的計(jì)算結(jié)果在Re較小時(shí)與實(shí)驗(yàn)值符合較好,而在Re較大時(shí)與實(shí)驗(yàn)值的偏差隨著Re的增大而變大;Cheng關(guān)系式的計(jì)算結(jié)果整體與實(shí)驗(yàn)值偏差較大,尤其是在過渡工況的起始Re附近偏差最大。因此,對(duì)于CEFR燃料組件這種幾何參數(shù)下的帶繞絲棒束組件,利用Engel關(guān)系式來進(jìn)行低Re下組件的阻力系數(shù)計(jì)算準(zhǔn)確程度較高。

圖3示出Re=534時(shí)與入口距離為0、1/3螺距、2/3螺距、1個(gè)螺距平面位置的速度場(chǎng)分布。可看出,由于采用平移周期性邊界條件,出入口的流場(chǎng)保持一致,可模擬充分發(fā)展段帶繞絲棒束組件的流動(dòng)。同時(shí),1/3螺距、2/3螺距平面位置的速度場(chǎng)分布可通過將入口速度場(chǎng)分別順時(shí)針旋轉(zhuǎn)120°、240°得到,這主要是繞絲順時(shí)針旋轉(zhuǎn)纏繞在燃料棒產(chǎn)生的效果。值得注意的是,由于燃料棒與燃料盒之間的間距較大使得流動(dòng)阻力較小,造成這些區(qū)域的冷卻劑流速最大,位于中心的燃料棒周圍的流速反而不大,這對(duì)于中心位置的燃料棒的冷卻是不利的。

2) 19根帶繞絲棒束組件的水力特性

19根帶繞絲棒束組件計(jì)算結(jié)果如圖4、5所示。與61根帶繞絲棒束組件的結(jié)果相同,利用SST轉(zhuǎn)捩模型計(jì)算得到的結(jié)果大多情況下與Engel關(guān)系式符合較好(差別小于5%),與ICRT關(guān)系式結(jié)果的差別稍大。而層流模型計(jì)算結(jié)果與經(jīng)驗(yàn)關(guān)系式差別較大。因此利用SST轉(zhuǎn)捩模型分析繞絲幾何參數(shù)一致的19根燃料棒組件的阻力系數(shù)也將與實(shí)際情況符合較好。

3) 7根帶繞絲棒束組件的水力特性

7根帶繞絲棒束組件計(jì)算結(jié)果如圖6、7所示。將層流模型和SST轉(zhuǎn)捩模型計(jì)算的結(jié)果與Engel關(guān)系式及ICRT模型結(jié)果比較發(fā)現(xiàn),SST轉(zhuǎn)捩模型計(jì)算的結(jié)果在大多情況下與Engel關(guān)系式符合較好(差別小于5%),僅在Re<300之后差別逐漸增大。而該計(jì)算結(jié)果與ICRT模型結(jié)果僅在Re=400附近較接近,隨著Re的增大或減小,這種差別逐漸增大。而利用層流模型計(jì)算的結(jié)果則與上述兩個(gè)關(guān)系式相差較大:層流模型計(jì)算結(jié)果與上述Engel關(guān)系式結(jié)果的差別保持在36%~48%,與ICRT關(guān)系式結(jié)果的差別保持在30%~36%。因此對(duì)于7根帶繞絲棒束組件的阻力系數(shù),SST轉(zhuǎn)捩模型的計(jì)算結(jié)果也與實(shí)際情況符合較好。

圖3 61根帶繞絲棒束組件不同截面速度分布Fig.3 Velocity distribution for different cross sections of 61 wire-wrapped rod bundle assembly

圖4 19根帶繞絲棒束組件水力特性計(jì)算結(jié)果與各關(guān)系式比較Fig.4 Comparison of hydraulic characteristic calculation result of 19 wire-wrapped rod bundle assembly with different relational expressions

綜上所述,通過將61、19及7根燃料棒的帶繞絲棒束組件的水力特性計(jì)算結(jié)果與實(shí)驗(yàn)以及經(jīng)驗(yàn)關(guān)系式對(duì)比發(fā)現(xiàn),利用SST轉(zhuǎn)捩模型計(jì)算低Re下組件水力特性的結(jié)果與實(shí)驗(yàn)及經(jīng)驗(yàn)關(guān)系式符合較好,因此該模型可應(yīng)用于低Re下此類繞絲結(jié)構(gòu)的分析計(jì)算。在后續(xù)的分析計(jì)算中均使用該模型。

4.2 4個(gè)螺距帶繞絲棒束組件

帶繞絲棒束組件的水力特性對(duì)于組件設(shè)計(jì)、堆芯流量分配以及事故工況下的流動(dòng)分析都具有重要意義。由于繞絲帶來的橫向流動(dòng),帶繞絲棒束組件的水力特性實(shí)驗(yàn)測(cè)量時(shí)需考慮橫向流動(dòng)帶來的影響。

由于4個(gè)螺距的帶繞絲棒束組件水力特性分析計(jì)算需要相當(dāng)龐大的網(wǎng)格數(shù)量,因此本文選用7根燃料棒的帶繞絲棒束組件進(jìn)行組件出、入口段及充分發(fā)展段水力特性的分析。

在4個(gè)螺距的7根帶繞絲棒束組件水力特性的計(jì)算中,計(jì)算模型采用SST轉(zhuǎn)捩模型。將組件盒壁面的6個(gè)面編號(hào),分別得到每個(gè)面與距離組件入口0、1/4螺距、1/3螺距、1/2螺距、2/3螺距、1個(gè)螺距、2個(gè)螺距、3個(gè)螺距、10/3螺距、11/3螺距、4個(gè)螺距的平面的交線中點(diǎn)的靜壓。組件各面編號(hào)及長(zhǎng)度方向位置(軸向各位置分別與圖中a~k字母對(duì)應(yīng))如圖8所示。

圖5 19根帶繞絲棒束組件不同截面速度分布Fig.5 Velocity distribution for different cross sections of 19 wire-wrapped rod bundle assembly

圖6 7根帶繞絲棒束組件水力特性計(jì)算結(jié)果與各關(guān)系式比較Fig.6 Comparison of hydraulic characteristic calculation result of 7 wire-wrapped rod bundle assembly with different relational expressions

1) 入口段壓力分布

比較組件6個(gè)組件盒壁面在a、c、e、f處(分別對(duì)應(yīng)0、1/3螺距、2/3螺距、1個(gè)螺距的位置)的靜壓,如圖9所示。從圖中可看出,在同一軸向位置,組件6個(gè)面上的壓力不一致,這也是由于繞絲帶來的橫向流動(dòng)產(chǎn)生的影響。這種各面壓力不一致的差別會(huì)隨著Re的增大而越發(fā)明顯[4]。位于a和f處各面的壓力分布趨勢(shì)基本一致,而兩個(gè)中間位置c和e的壓力分布有所不同。

實(shí)際上,由于繞絲周期性的纏繞,各面的壓力分布也存在著周期性的旋轉(zhuǎn)變化。對(duì)圖9中

圖7 7根帶繞絲棒束組件不同截面速度分布Fig.7 Velocity distribution for different cross sections of 7 wire-wrapped rod bundle assembly

數(shù)據(jù)稍加處理,即將c、e面上6個(gè)點(diǎn)與a或f面上6個(gè)點(diǎn)根據(jù)繞絲的旋轉(zhuǎn)方向進(jìn)行對(duì)應(yīng)(類似“移動(dòng)相位”),即c面上點(diǎn)3對(duì)應(yīng)a面上的點(diǎn)1,e面上點(diǎn)5對(duì)應(yīng)a面上點(diǎn)1,以此類推。同時(shí)將每個(gè)面上各點(diǎn)的壓力換算成與該面各點(diǎn)壓力平均值的差(p-p平面平均),處理后的結(jié)果如圖10所示。從圖10可看出,經(jīng)調(diào)整后各壁面的壓力的變化趨勢(shì)一致,即由于周期性橫向流動(dòng)的影響,不同軸向位置的壓力分布與速度分布同樣是旋轉(zhuǎn)對(duì)應(yīng)的。注意到e和f平面的壓力分布曲線進(jìn)行相位移動(dòng)后基本重合,表明e、f位置處的流動(dòng)已進(jìn)入充分發(fā)展階段。

圖8 組件壁面和軸向位置Fig.8 Position on wall of assembly and vertical direction

圖9 入口段組件各壁面壓力分布Fig.9 Pressure distribution on different walls at inlet of assembly

2) 出口段壓力分布

出口段壓力分布與入口段的類似,同樣列出各點(diǎn)的壓力分布和進(jìn)行調(diào)整處理后的分布如圖11、12所示。從圖12可看出,經(jīng)調(diào)整后各壁面的壓力的變化趨勢(shì)一致,即由于周期性橫向流動(dòng)的影響,不同軸向位置的壓力分布與速度分布同樣是旋轉(zhuǎn)對(duì)應(yīng)的。注意到各平面的壓力分布曲線進(jìn)行相位移動(dòng)后基本重合,表明出口段的影響很小。

圖10 調(diào)整后的入口段組件各壁面壓力分布Fig.10 Pressure distribution on different wallsat inlet of assembly after phase change

圖11 出口段組件各壁面壓力分布Fig.11 Pressure distribution on different walls at outlet of assembly

圖12 調(diào)整后的出口段組件各壁面壓力分布Fig.12 Pressure distribution on different wallsat outlet of assembly after phase change

3) 整數(shù)倍螺距處壓力分布

圖13為組件上各整數(shù)倍螺距處的壓力分布。從圖13可發(fā)現(xiàn),從第1個(gè)螺距位置的壓力分布開始到第4個(gè)螺距位置的壓力分布曲線相互平行,意味著流動(dòng)進(jìn)入到充分發(fā)展。相互平行的曲線意味著同一面上1個(gè)螺距的壓降均相等,而且與其他面上1個(gè)的螺距的壓降一致,并且同一面上整數(shù)倍螺距的壓降恰好是1個(gè)螺距壓降的相應(yīng)整數(shù)倍,這說明雖然繞絲產(chǎn)生的橫向流動(dòng)使組件6個(gè)壁面上壓力分布有所不同,但同一壁面上壓降沿著軸向按螺距均勻分布。這樣在充分發(fā)展段,組件的壓降沿著軸向也是按螺距均勻分布的,組件每個(gè)螺距的阻力系數(shù)不變。這一結(jié)論也與文獻(xiàn)[16]中的實(shí)驗(yàn)結(jié)果相符。在進(jìn)行帶繞絲棒束組件水力特性測(cè)量時(shí),需在組件同一面上按照整數(shù)倍螺距來布置測(cè)點(diǎn),才能避免由于橫向流動(dòng)對(duì)測(cè)量帶來的影響。

圖13 整數(shù)倍螺距處組件各壁面壓力分布Fig.13 Pressure distribution on different walls at positions of integer multiples of helix

通過上述結(jié)果可知,在較低Re下該帶繞絲棒束組件的入口穩(wěn)定段長(zhǎng)度小于1/2螺距的,出口段的影響很小,其余部分為充分發(fā)展段。在充分發(fā)展段內(nèi)部,流體流動(dòng)按照螺距呈周期性分布,因此,在分析計(jì)算該類帶繞絲棒束組件時(shí),只需保證入口、出口段以及至少1個(gè)螺距的充分發(fā)展段即可充分說明組件內(nèi)流體的流動(dòng)情況,其余為周期性部分。因此本文中4個(gè)螺距的計(jì)算結(jié)果能涵蓋全尺寸組件。

5 結(jié)論

1) 通過將CFX計(jì)算結(jié)果與俄羅斯實(shí)驗(yàn)數(shù)據(jù)以及國(guó)際上的相關(guān)經(jīng)驗(yàn)關(guān)系式比較可得到,利用SST轉(zhuǎn)捩模型計(jì)算帶繞絲棒束組件低Re下的水力特性可得到比較好的結(jié)果,因此利用該模型計(jì)算低Re下CEFR燃料組件的水力特性是合適的。

2) 通過對(duì)4個(gè)螺距的7根帶繞絲棒束組件的水力特性分析表明,在較低Re下該帶繞絲棒束組件的入口穩(wěn)定段長(zhǎng)度小于1/2螺距的,出口段的影響很小。同時(shí),在流動(dòng)達(dá)到充分發(fā)展后,雖然繞絲產(chǎn)生的橫向流動(dòng)使組件6個(gè)壁面上壓力分布有所不同,但同一壁面上壓降沿著軸向按螺距是均勻分布的,從而組件的壓降沿著軸向也是按螺距均勻分布的,組件每個(gè)螺距的阻力系數(shù)不變。

3) 計(jì)算結(jié)果給實(shí)驗(yàn)的指導(dǎo):在進(jìn)行帶繞絲棒束組件水力特性測(cè)量時(shí),需在組件同一面上按照整數(shù)倍螺距來布置測(cè)點(diǎn),才能避免由于橫向流動(dòng)對(duì)測(cè)量帶來的影響。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 人妻一本久道久久综合久久鬼色| 成年人福利视频| 日韩精品无码免费一区二区三区| 欧美影院久久| 亚洲国产日韩一区| 欧美激情第一欧美在线| 国产美女在线观看| 台湾AV国片精品女同性| 亚洲人成高清| av天堂最新版在线| 久久久久中文字幕精品视频| 国产h视频在线观看视频| 国产国语一级毛片在线视频| 啦啦啦网站在线观看a毛片 | 亚洲bt欧美bt精品| 国产女人18水真多毛片18精品 | 成人免费午夜视频| 国产免费看久久久| 日韩欧美国产成人| 东京热一区二区三区无码视频| 欧洲成人在线观看| 国产成人久久777777| 手机在线免费不卡一区二| 国产乱视频网站| 永久免费无码日韩视频| 免费在线色| 成人一区专区在线观看| 欧美在线观看不卡| 国产精品综合久久久| 成人福利在线观看| 伊人网址在线| 青青草原国产av福利网站| 成人91在线| 国产精品性| 欧美亚洲国产精品第一页| 91无码网站| 亚洲欧美另类日本| 波多野结衣中文字幕久久| 91麻豆国产在线| 天堂成人在线视频| 日韩久草视频| 青青青视频91在线 | 一级全免费视频播放| 国产亚洲精久久久久久久91| 亚洲国产综合精品中文第一| 亚洲色无码专线精品观看| 波多野结衣第一页| 91视频青青草| 91精品久久久久久无码人妻| 四虎成人精品| 免费毛片全部不收费的| 日本午夜网站| 成人在线综合| 亚洲大尺码专区影院| 国产久操视频| 中文无码伦av中文字幕| 国产乱子伦精品视频| 精品成人一区二区三区电影| 久久久久久久蜜桃| 国产日韩欧美精品区性色| 国产午夜精品一区二区三区软件| 国产三区二区| 亚洲综合精品香蕉久久网| 欧美一级一级做性视频| 欧洲在线免费视频| 黄片在线永久| 国产一区二区人大臿蕉香蕉| 国产精品.com| 欧美精品啪啪一区二区三区| 欧美另类精品一区二区三区| 无码网站免费观看| 欧美日韩国产在线播放| 国产二级毛片| 国产人人乐人人爱| 2021天堂在线亚洲精品专区 | 欧洲亚洲一区| 国产一区二区三区在线精品专区| 大香伊人久久| 男女精品视频| 91精品网站| 亚洲色欲色欲www网| 国产乱人伦偷精品视频AAA|