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打印中的模型分割與打包
主站蜘蛛池模板: 毛片三级在线观看| 国产成人一区免费观看| 国产小视频免费| 欧美日韩中文国产va另类| 亚洲日韩国产精品无码专区| 国产乱肥老妇精品视频| 亚洲九九视频| 九九精品在线观看| 91色在线视频| 免费人成在线观看成人片 | 久久公开视频| 狠狠色噜噜狠狠狠狠色综合久| 亚洲AV人人澡人人双人| 波多野结衣一区二区三区四区 | 日韩在线网址| 日本爱爱精品一区二区| 538国产视频| 亚洲精选高清无码| 亚洲黄网视频| 婷婷色一区二区三区| 最新国产麻豆aⅴ精品无| 好紧好深好大乳无码中文字幕| 中文字幕调教一区二区视频| 日韩国产一区二区三区无码| 美女潮喷出白浆在线观看视频| 精品国产免费观看一区| 久久久受www免费人成| 国产精品专区第一页在线观看| 亚洲天堂网在线播放| 国产精品妖精视频| 一本视频精品中文字幕| 免费A级毛片无码免费视频| 亚洲av无码人妻| 免费午夜无码18禁无码影院| 五月激激激综合网色播免费| 欧洲亚洲欧美国产日本高清| 亚洲欧美色中文字幕| 欧美a级在线| 毛片最新网址| 亚洲精品制服丝袜二区| 日韩av无码精品专区| 国产女人在线观看| 久久女人网| a级毛片网| 99re在线观看视频| 亚洲视频无码| 国产国语一级毛片在线视频| 久久人体视频| 亚洲一级毛片在线观播放| 女同久久精品国产99国| 日韩免费中文字幕| 91国内视频在线观看| 亚洲国产中文精品va在线播放| 国产精品3p视频| 青青青草国产| 国产丝袜91| 成年女人a毛片免费视频| 黄片在线永久| 乱系列中文字幕在线视频| 五月婷婷激情四射| 国产真实乱了在线播放| www.av男人.com| 亚洲一级毛片免费观看| 伊人成人在线| 日本一本在线视频| 丁香亚洲综合五月天婷婷| 国产精品吹潮在线观看中文| 欧美精品H在线播放| 日韩av在线直播| 欧美精品H在线播放| 毛片网站在线播放| 欧美a级完整在线观看| 久久永久精品免费视频| 在线观看欧美国产| 一级高清毛片免费a级高清毛片| 91热爆在线| 91麻豆精品国产高清在线| 日本成人在线不卡视频| 丁香六月激情婷婷| 久久香蕉国产线| 国产精品乱偷免费视频| 亚洲精品午夜无码电影网|