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

基于縮比實驗模型的快堆主容器內氬氣空間傳熱特性數值研究

2023-12-16 05:44:16羅心越陸道綱馮佳琪于宗玉張鈺浩
核科學與工程 2023年5期
關鍵詞:實驗

羅心越,陸道綱,馮佳琪,于宗玉,張鈺浩,*

基于縮比實驗模型的快堆主容器內氬氣空間傳熱特性數值研究

羅心越1,2,陸道綱1,2,馮佳琪1,2,于宗玉1,2,張鈺浩1,2,*

(1. 華北電力大學核科學與工程學院,北京 102206;2. 非能動核能安全技術北京市重點實驗室,北京 102206)

池式鈉冷快堆主容器內鈉液面上方覆蓋有氬氣,鈉液面通過多種傳熱方式相互耦合,不斷向堆容器上部結構傳遞熱量。為了對上部結構進行溫度載荷下的強度分析,有必要開展氬氣空間傳熱特性研究。一方面,國內外針對錐頂蓋氬氣空間傳熱特性開展了較多的數值模擬研究,但大多缺乏實驗驗證。另一方面,雖有少量針對氬氣空間的換熱特性實驗研究,但由于測量技術的制約,僅能獲得關鍵部件的溫度分布,難以獲得其流動特性,因此對氬氣空間復雜的傳熱特性和傳熱機理有待進一步研究。本文基于快堆鈉液面和氬氣空間換熱特性縮比實驗,開展氬氣空間熱工水力特性數值模擬,通過與實驗結果對比,驗證了數值方法的有效性,并在此基礎上,對其溫度分布及流動特性進行詳細研究。研究結果表明:在不同工況下,氬氣空間內溫度分布總體對稱,由于徑向熱屏蔽的存在,其上部窄流道區域溫度明顯低于下部氬氣大空間溫度;徑向熱屏蔽與底部鈉液面間的氬氣大空間內,由于瑞利()數較高,自然對流現象強烈,換熱效果較強,在靠近底部鈉液面處形成局部渦流,而徑向熱屏蔽外部窄流道區域內由于數較小,自然對流較弱,主要依靠熱輻射和導熱進行傳熱。

鈉冷快堆;錐頂蓋;氬氣空間;對流換熱

池式鈉冷快堆是第四代核電技術中的重要堆型,其堆本體結構復雜,包含旋塞、氬氣空間、中間熱交換器(IHX)、獨立熱交換器(DHX)、泵支承、溫度液位計等多個部件。堆本體內一回路采用液態鈉作為冷卻劑,錐頂蓋空間內覆蓋氬氣作為保護氣體。反應堆運行時,鈉液面不斷向堆容器的上部結構(徑向熱屏蔽、斜肩等)及氬氣傳遞熱量。錐頂蓋氬氣空間內傳熱機理復雜,存在多種傳熱方式相互耦合的復雜行為,錐頂蓋溫度場分布難以準確獲得。此外,反應堆運行時,液鈉蒸氣在氬氣空間內形成的氣溶膠對熱輻射具有吸收作用,錐頂蓋、旋轉屏蔽塞和貫穿件支承在徑向和軸向方向上存在局部熱傳導等,這進一步增加了錐頂蓋溫度場分布的復雜性,給錐頂蓋溫度載荷下的強度分析和應力評價及設計帶來困難。

針對氬氣空間換熱特性,部分學者開展了基于簡化模型的實驗研究。吳強等[1]建立模擬快堆氬氣空間的水-空氣實驗裝置,研究封閉空間內氣體溫度分布特征,發現液面上存在兩種換熱方式-輻射換熱和對流換熱,且液面的蒸發現象能夠阻礙兩種換熱效果。李天舒等[2]基于相似原理建立快堆氬氣空間的簡化實驗裝置,實驗結果顯示氣體溫度在封閉空間中間大部分區域維持相對恒定,同時提出了氬氣空間高度和主容器冷卻系統設計應統籌考慮的設計思想。

為了探究錐頂蓋及氬氣空間復雜熱工水力特性,國內一些學者建立了詳細的三維模型,張曉露[7]針對錐頂蓋主要部件結構進行精細建模,分析旋塞與旋塞支承頸氣隙及生物屏蔽柱對金屬層溫度分布的影響,發現生物屏蔽柱對金屬層的輻射效果較強,30 mm氣隙對錐頂蓋整體的保溫效果最好。余新太等[8]應用低雷諾數湍流模型和Discrete Ordinates(DO)輻射傳熱模型對正常運行工況下示范快堆氬氣空間進行了數值模擬,計算結果表明,泵支承、中間熱交換器、獨立熱交換器周向溫差較大,需注意彎曲變形。林超等[9]對氬氣空間及相關結構部件進行穩態共軛傳熱計算,發現原型設計方案泵支承結構位于氬氣空間部分周向溫度分布有較大溫度梯度,泵支承結構靠近主容器支承頸側溫度明顯高于遠離支承頸側。盡管以上數值模擬揭示了氬氣空間內部部分復雜傳熱行為和機理特性,但其結果大多缺乏實驗驗證,針對數值模擬中的參數選取、計算方法也存在較大不確定性。因此,開展主容器頂部空間內氬氣的復雜流動特性的實驗研究十分必要。

綜上,由于池式快堆頂部氬氣空間結構復雜、耦合傳熱機理復雜且受到多種外部因素的影響,難以開展等比例實驗研究主容器頂部空間內氬氣的復雜流動特性,因此,為了揭示氬氣空間的換熱機理、驗證其換熱模型,基于相似原理建立了快堆鈉液面和氬氣空間換熱縮比實驗臺架,同時,鈉液面具有流動性,其加熱所產生的鈉蒸氣會對氬氣空間內部氬氣的熱工水力特性產生影響,為了排除鈉蒸氣氣溶膠對輻射換熱效果的影響,開展了采用不銹鋼模擬鈉液面運行工況及換料工況下的傳熱特性研究實驗,并開展池式快堆主容器頂部氬氣空間熱工水力特性數值模擬,驗證計算方法的有效性,進而探究氬氣空間內多種傳熱方式并存的復雜耦合傳熱機理,為主容器頂部空間溫度載荷下的強度分析應力評定提供關鍵參數。

1 實驗簡介

在快堆主容器內氬氣空間換熱理論分析基礎上,搭建鈉液面與氬氣空間換熱縮比實驗臺架。實驗裝置主要包括通風系統、電加熱系統、實驗裝置以及數據采集系統。其中,實驗裝置基本保留了氬氣空間的主要結構,包含錐頂蓋、獨立熱交換器(DHX)、中間熱交換器(IHX)、泵支承、徑向熱屏蔽、旋塞等。底部為封閉的不銹鋼圓筒,加熱裝置在不銹鋼底板的底面以完成實驗裝置的加熱,實驗裝置結構如圖1所示。因此該實驗裝置可用于模擬有/無鈉液面條件下快堆主容器頂部空間覆蓋氣體的自然循環特性。

圖1 實驗臺架整體結構圖

測量及數據采集系統,包括熱電偶、壓力表、流量計及數據采集系統(NI PIXe),用來進行數據采集、存儲和后處理。其中,通過將熱電偶點焊在錐頂蓋表面以測量其表面溫度,將熱電偶布置在不同高度的支架上以監測氬氣空間內溫度分布,其位置如圖2所示。每個支架上隨軸向高度布置6個測點。

圖2 測溫點分布示意圖

2 數值模型及方法

2.1 幾何模型

基于鈉液面與氬氣空間換熱縮比實驗臺架建立數值計算幾何模型,如圖3所示。

圖3 氬氣空間三維模型結構

本研究使用STAR-CCM+軟件對主容器錐頂蓋進行建模。錐頂蓋內包含獨立熱交換器(DHX)、中間熱交換器(IHX)、泵支承、徑向熱屏蔽、旋塞等多個部件,同時存在許多窄縫流道。為了能夠實現模型關鍵部件的主要運行功能,模擬鈉冷快堆頂部氬氣空間的自然循環現象,同時考慮到整體網格質量、生成難度以及保證計算精度,根據結構特點和流動特性,對各部件單獨進行網格劃分,將氬氣空間分為貫穿件、錐頂蓋金屬層保溫層、旋塞、凸臺等不同模塊。同時,在建模過程中盡可能保留錐頂蓋的主要幾何特征,將螺栓、熱電偶、液位計等尺寸較小的部件省略,將貫穿件及其周圍屏蔽簡化為實體或空心圓柱,將旋塞屏蔽層簡化為實體結構,以便有效模擬錐頂蓋氬氣空間內的關鍵熱工特性。

基于所建立的幾何模型開展非結構化網格劃分,并進行網格敏感性分析,由于鈉冷快堆頂部氬氣空間內部結構復雜,建立2 500萬、4 900萬、8 000萬的網格開展網格無關性分析計算,對比錐頂蓋平均溫度,結果如圖4所示。

圖4 網格無關性分析

由圖4可見,2 500萬網格和4 900萬網格模擬結果相差較大,4 900萬網格與8 000萬網格的模擬結果差異小于3%,因此,能夠認為,4 900萬網格計算結果滿足計算要求。最終選取約4 900萬網格開展驗證分析。

2.2 數值模型及邊界條件

2.2.1數值模型

在近壁面區域,湍流應力不起作用,因此,在近壁面區域,采用壁面函數法,無需將壁面區域加密,用一組半經驗公式將壁面上的物理量與湍流核心區內相應物理量相聯系,直接得到與壁面相鄰控制體積的節點變量值[11],在本文中,采用STAR-CCM+默認的壁面模型Ally+wall treatment計算。

DO模型能夠較好地解決絕大多數光學深度區間的輻射問題,能夠考慮到氣體介質與顆粒間的輻射交換,能量不平衡率最小,在存在局部熱源的情況下,能夠較好地反應模型輻射換熱情況。

本文選擇DO模型以求解有限數量的離散立體角的輻射運輸方程,其原理是通過計算每個離散縱坐標上的偏微分方程來求解強度[9],輻射傳熱方程形式如下:

——射線路徑長度;

——當地溫度。

DO模型是將球坐標系下的立體角轉化為直角坐標系下的方向,把沿方向的輻射傳熱方程視為某個場方程,于是式(3)變形為:

2.2.2計算工況及邊界條件

為了驗證不同工況下主容器頂部空間的傳熱特性,在實驗臺架中,針對不銹鋼底面模擬鈉液面溫度開展運行工況及換料工況兩組關鍵工況的實驗,在模擬計算中,不同工況下模擬鈉液面溫度的設置與實驗臺架保持一致。分別為:

(1)運行溫度工況:不銹鋼底面(模擬鈉液面)溫度為540 ℃,通過分析該工況下各堆內構件的關鍵三維熱動水力參數,評估在反應堆正常運行下頂部空間的壓力邊界安全性。

(2)換料溫度工況:不銹鋼底面(模擬鈉液面)溫度為250 ℃,基于該工況下氬氣空間的溫度分布特性,驗證在裝卸料情況下,反應堆頂部空間的壓力邊界完整性。

不同工況下邊界條件如表1所示。

表1 各位置邊界條件

基于上述關鍵邊界條件的設置,開展鈉液面與氬氣空間換熱實驗數值模擬計算。計算過程中所涉及的不銹鋼、保溫層及氬氣的密度、導熱系數、比熱容等物性參數均來源于《載熱質熱物性計算程序及數據手冊》[10]。

3 氬氣空間傳熱特性分析

3.1 數值方法驗證

在氬氣空間數值模擬中,運行工況、換料工況下錐頂蓋溫度分布如圖5所示。

圖5 錐頂蓋溫度分布

如圖5所示,錐頂蓋上溫度總體呈對稱分布,運行工況下,錐頂蓋溫度基本穩定在260~400 ℃,換料工況下,錐頂蓋溫度穩定在80~120 ℃。整體溫度分布較均勻,但在貫穿件附近,溫度出現明顯上升,推測因為氬氣與貫穿件的對流換熱加強,導致靠近貫穿件側的部分溫度升高。

對照實驗臺架布置的測點位置,在數值模擬計算中選擇相同位置處(見圖3)的溫度進行對比分析。選擇Line 2作為監測線,對比如圖6所示,分別分析實驗運行工況、模擬運行工況、實驗換料工況、模擬換料工況的溫度變化。

圖6 Line 2溫度變化曲線

通過溫度變化對比曲線可以看出,數值模擬計算結果的溫度變化趨勢與實驗結果基本一致,與實驗結果相比,數值模擬得到溫度平均相對誤差在 5%以內。可以驗證,數值模擬所采用的計算模型、計算方法合理,能夠較好地描述主容器頂部氬氣空間內覆蓋氣體的熱工水力特性。

3.2 氬氣空間傳熱特性分析

3.2.1氬氣空間溫度分布

在氬氣空間內選擇監測面截面(1)、截面(2),得到氬氣空間的溫度分布如圖7所示。可以發現,氬氣空間的溫度分布基本呈對稱狀態,運行工況下最高溫度為 540 ℃,最低溫度約為411 ℃,溫差約 131 ℃,換料工況最高溫度為250 ℃,最低溫度約為 151 ℃,溫差約 99 ℃,垂直高度上溫度自下而上逐漸降低,最高溫度出現在鈉液面附近。由于徑向熱屏蔽的存在,上部窄流道區域溫度明顯低于下部氬氣大空間溫度,其存在非常有效的阻隔了溫度的傳導,使金屬層和保溫層的溫度始終處于主容器的設計限值內。

圖7 監測面溫度分布

此外,從截面溫度云圖中能夠明顯看出,在氬氣空間底部存在兩個對稱位置的局部渦流區域,在該區域內,氬氣溫度分布呈明顯的不規律性,形成兩個小的漩渦。提取換料溫度工況下截面速度矢量分布如圖8所示。靠近通道屏蔽的一側渦流區域,由于氬氣在該位置同時與通道屏蔽和金屬層進行對流傳熱和輻射換熱,溫度明顯較低,密度增大,向下回流,從而形成局部渦流。

圖8 氬氣空間整體流場

在數值模擬中選擇與實驗裝置相同位置處(見圖 3)的監測線進行溫度對比分析,如圖9所示。

基于線上各點的溫度分布可知,氬氣空間內溫度隨垂直高度的增加而降低,變化基本趨勢不變,Line 2、Line 3在氬氣空間內呈對稱狀態,其溫度變化趨勢基本一致,進一步驗證在整個氬氣空間內,溫度分布基本呈對稱狀態,分析Line 4可以看到,在貫穿件附近,由于通道屏蔽和凸臺的存在,氬氣的自然循環高度降低,因此,相較Line 1,Line 4的變化趨勢更平緩、溫度梯度更小。而由于氬氣與貫穿件的熱交換作用明顯,氬氣溫度會出現小范圍波動。

3.2.2換熱機理分析

在傳熱過程中,氬氣在空間內流動,將熱量傳遞給徑向熱屏蔽,經熱傳導傳遞到徑向熱屏蔽外表面,最后通過外表面傳遞到金屬層。自然循環過程中,自然對流傳熱是氬氣空間內引起氬氣與物體表面之間的熱量傳遞方式。

圖9 不同工況下監測線溫度曲線圖

瑞利()數是與自然對流相關的無量綱數,反映自然對流效果,其公式為:

式中,格拉曉夫()數描述了流體的浮力和粘度之間的關系,普朗特數()描述了動量擴散系數和熱擴散系數之間的關系。普朗特數()反映了動量擴散厚度和能量擴散厚度的相對大小。

——重力加速度為9.8 m/s2;

為了分析氬氣空間內基于自然對流機理下的氬氣換熱特性,在氬氣空間內選擇垂直于徑向熱屏蔽的兩條分段平行的溫度線,如圖2所示。提取各監測線溫度,計算對應瑞利數如表2所示。

表2 瑞利(Ra) 數表

圖10為不同工況下溫度線在無量綱高度上的溫度分布。從圖中可知,溫度隨無量綱高度增加逐漸降低,無量綱高度為0~0.2,溫度變化劇烈且變化梯度大,無量綱高度在0.2~1,溫度下降趨勢變緩,結合瑞利()數分析可得,在徑向熱屏蔽下氬氣空間大區域內,數達到105~107量級,空間越大,瑞利()數越大,自然對流特征越明顯,驅動力越大,因此會產生局部渦旋(見圖7)。而在徑向熱屏蔽外部窄流道區域內,數較小,自然對流影響較小,溫度分布更均勻,主要依靠輻射和導熱進行換熱。

圖10 監測位置溫度分布變化

4 結論

本文基于快堆鈉液面和氬氣空間換熱臺架運行工況及換料工況下的傳熱特性研究實驗,對鈉冷快堆頂部氬氣空間開展整體三維數值模擬分析。獲得主要結論如下:

(2)氬氣空間溫度分布整體基本呈對稱狀態,在垂直高度上自下而上出現明顯溫度分層,平均溫差達到150 ℃。由于徑向熱屏蔽有效阻隔了熱量的傳導,上部窄流道區域溫度明顯低于下部氬氣大空間溫度。

(3)在徑向熱屏蔽下的氬氣空間大區域內,數能夠達到105~107量級,自然對流換熱特征明顯,形成局部渦流,但在徑向熱屏蔽外部窄流道區域,數較小,自然對流換熱影響較小,主要依靠輻射和導熱進行換熱。

[1] 吳強,張東輝,吳純良,等. CEFR主容器氬氣空間溫度分布模擬實驗研究[G]. 中國原子能科學研究院年報,2009:9-10.

[2] 吳強,張東輝,李天舒,等. 基于相似理論的快堆氬氣空間溫度場實驗研究[J]. 核科學與工程,2016,36 (04):465-469.

[3] 吳強,張東輝,劉云焰. 快堆主容器內氬氣空間傳熱實驗原理及分析[J]. 核科學與工程,2009,29(01):27-32.

[4] 許義軍. CEFR大小旋塞氬氣空間流固耦合傳熱計算的數值分析[C]. 計算流體力學研究進展——第十二屆全國計算流體力學會議論文集,2004:726-731.

[5] 張維忠,賈斗南,秋穗正. 高溫氣冷堆一回路艙室的輻射換熱計算模型[J]. 核動力工程,2002(02):10-16.

[6] 郭磊. 對FLUENT輻射模型的數值計算與分析[J]. 制冷與空調(四川),2014,28(03):358-360.

[7] 張曉露. 快堆主容器頂部復雜空間中的覆蓋氣體熱工水力行為研究[D]. 北京:華北電力大學(北京),2019.

[8] 余新太,高鑫釗,馬曉,等. 示范快堆主容器內氬氣空間數值模擬[J]. 核科學與工程,2021,41(04):695-702.

[9] 林超,楊紅義,周志偉. CFR600泵支承氬氣空間部分三維傳熱數值模擬[J]. 原子能科學技術,2021,55(04):654-659.

[10]居懷明,徐元輝,李懷萱. 載熱質熱物性計算程序及數據手冊[M]. 北京:原子能出版社,1990.

[11]謝東. 地下水電站廠房氣流組織CFD數值模擬方法研究[D]. 重慶:重慶大學,2015.

[12]李光正,馬洪林. 封閉腔內高瑞利數層流自然對流數值模擬[J]. 華中科技大學學報(城市科學版),2004(03):14-17.

[13]邱志民. AP1000非能動余熱排出系統一、二次側耦合傳熱特性實驗研究[D]. 北京:華北電力大學(北京),2021.

[14]梁江濤,陸道綱,趙海琦,等. 池式鈉冷快堆雙環路12%差異非對稱功率運行及流量調節緩解工況的三維數值模擬[J]. 科學技術與工程,2021,21(25):10703-10709.

[15] Dinesh Nath,Mahendra K.Verma.Numerical simulation of convection of argon gas in fast breeder reactor[G]. Annals of Nuclear Energy,2014:63.

[16] Xiaoxue Huang,Shuisheng He,Numerical modelling of cover gas thermal hydraulics in Sodium-cooled Fast Reactors[J]. Nuclear Engineering and Design,2019,355:110347.

[17] Dinesh Nath,Mahendra K.Verma.Numerical simulation of convection of argon gas in fast breeder reactor[G]. Annals of Nuclear Energy,2014:63.

[18] Guleren K M,Pinarbasi A. Numerical simulation of the stalled flow within a vaned centrifugal pump. Proceedings of the Institution of Mechanical Engineers,Part C[J]. Journal of Mechanical Engineering Science,2004;218(4):425-435.

Numerical Simulation of the Argon Space for the Heat Transfer Characteristics in the Main Vessel Based on the Scaled-down Test Model

LUO Xinyue1,2,LU Daogang1,2,FENG Jiaqi1,2,YU Zongyu1,2,ZHANG Yuhao1,2,*

(1. School of Nuclear Science and Engineering,North China Electric Power University,Beijing 102206,China;2. Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy,Beijing 102206,China)

The argon gas is covered above the sodium level in the main vessel of pool-type sodium-cooled fast reactor (SFR). The sodium level continuously transfers heat to upper structure of the roof slab through the coupled heat transfer mechanism. In order to analyze the strength and evaluate the stress of roof slab, it is necessary to study the heat transfer characteristics of the argon space. On the one hand, many numerical simulation studies have been carried out on the space heat transfer characteristics of the argon gas in the roof slab. Most of the numerical simulation studies on the roof slab and the argon space lack of experimental verification. In addition, due to the limitation of measurement technique, it is easy to obtain the temperature characteristics, but it is difficult to obtain the flow characteristics in the existing study of numerical simulation. Therefore, the complex heat transfer characteristics and mechanism of the argon space still needs study. In this paper, a numerical simulation of the thermal-hydraulic characteristics in the argon space was carried out based on the heat-transfer experimental facility for the argon space at the top of the main vessel (HEFA), whose results were compared with the experimental data. The results show that the numerical method is effective, and the calculated results are reasonable. By analyzing the three-dimensional temperature distribution of the domain, it is found that under different working conditions, the temperature distribution in the argon space is generally symmetric. Due to the existence of radial heat shielding, the temperature in the upper narrow region is obviously lower than that in the external narrow channel area. In addition, there is a vortex forming in the local argon space. The heat transfer effect of natural convection is obvious, while thenumber of argon near the radial thermal shield is small in the external narrow channel area. The main heat transfer processes are heat conduction and radiation heat transfer.

Sodium cooled fast reactor; Roof slab; Argon space; Convective heat transfer

TL33

A

0258-0918(2023)05-0961-10

2022-10-04

羅心越(1996—),女,甘肅天水人,碩士研究生,現主要從事反應堆熱工水力學方面的研究

張鈺浩,E-mail:zhangyuhao@ncepu.edu.cn

猜你喜歡
實驗
我做了一項小實驗
記住“三個字”,寫好小實驗
我做了一項小實驗
我做了一項小實驗
記一次有趣的實驗
有趣的實驗
小主人報(2022年4期)2022-08-09 08:52:06
微型實驗里看“燃燒”
做個怪怪長實驗
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 亚洲男人在线天堂| 波多野结衣爽到高潮漏水大喷| 亚洲日本www| 天堂亚洲网| 制服丝袜在线视频香蕉| 国产xx在线观看| 国内精品手机在线观看视频| 国产视频资源在线观看| 亚洲天堂视频在线观看免费| 欧亚日韩Av| 人妻丰满熟妇av五码区| 国模粉嫩小泬视频在线观看 | 色窝窝免费一区二区三区| 亚洲欧美综合精品久久成人网| 美女国内精品自产拍在线播放| 欧美综合区自拍亚洲综合天堂 | 青青青视频免费一区二区| 午夜综合网| 99热这里只有精品国产99| 热思思久久免费视频| 亚洲高清无码久久久| 久久久亚洲国产美女国产盗摄| 精品成人免费自拍视频| 91欧美在线| 婷婷亚洲视频| 伊人色在线视频| 免费a级毛片视频| 91色国产在线| 欧美亚洲欧美| 尤物在线观看乱码| 久久久噜噜噜久久中文字幕色伊伊 | 亚洲日韩高清无码| 99er精品视频| 日韩久久精品无码aV| 欧美97欧美综合色伦图| 久久久久青草线综合超碰| 无码有码中文字幕| 呦女亚洲一区精品| 丁香五月亚洲综合在线 | 欧美日韩一区二区在线免费观看| 亚洲国产午夜精华无码福利| 日韩经典精品无码一区二区| 国产一二三区在线| 亚洲精品第1页| 亚洲综合18p| 亚洲天堂2014| 99激情网| 国产不卡在线看| a级高清毛片| 亚洲天堂日韩av电影| 婷婷伊人久久| 午夜日韩久久影院| 日本精品影院| 国产成年女人特黄特色毛片免 | 色哟哟国产精品| 99在线视频网站| 国产区成人精品视频| 无码专区国产精品第一页| 亚洲欧洲日韩久久狠狠爱 | 精品一区二区三区四区五区| 国产免费怡红院视频| 欧美在线伊人| 精品人妻一区二区三区蜜桃AⅤ| 伊人久久久久久久| 精品第一国产综合精品Aⅴ| 黄片一区二区三区| 国产精品伦视频观看免费| 五月婷婷亚洲综合| 欧美福利在线| 老司机午夜精品网站在线观看| 久视频免费精品6| 国产区福利小视频在线观看尤物| 午夜福利亚洲精品| 亚洲男女在线| 欧美另类图片视频无弹跳第一页| 欧美日韩第二页| 青青草原国产精品啪啪视频| 国产精品2| 亚洲中文字幕av无码区| 欧美日韩v| 久久免费视频6| 亚洲人成人伊人成综合网无码|