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

跨臨界泄壓瞬態(tài)傳熱特性的模型敏感性分析

2018-01-08 04:58:09宋美琪劉曉晶
核科學與工程 2017年6期

宋美琪,劉曉晶,程 旭

(上海交通大學核科學與工程學院,上海 200240)

跨臨界泄壓瞬態(tài)傳熱特性的模型敏感性分析

宋美琪,劉曉晶,程 旭

(上海交通大學核科學與工程學院,上海 200240)

在停堆或失水等事故工況下,超臨界水冷堆將經(jīng)歷跨臨界泄壓過程,系統(tǒng)壓力從超臨界狀態(tài)降到擬臨界點22.064MPa以下。而對于次臨界區(qū),臨界點附近的臨界熱流密度值很低,極易發(fā)生沸騰臨界,導(dǎo)致加熱棒壁面溫度迅速升高,因此跨臨界泄壓過程是超臨界水冷堆失水事故安全分析的關(guān)鍵。目前,跨臨界泄壓瞬態(tài)過程可以通過系統(tǒng)程序進行計算,但依然缺乏有效的實驗驗證。故本文依托上海交通大學的超臨界流體多功能實驗回路(Supercritical WAter MUltiPurpose loop,SWAMUP)跨臨界泄壓過程的實驗,利用德國核安全中心(GRS)開發(fā)的系統(tǒng)程序ATHLET3.0進行建模計算,分析跨臨界泄壓過程傳熱特性。通過調(diào)節(jié)次臨界區(qū)臨界熱流密度、最小膜態(tài)沸騰溫度、驟冷前沿模型等相關(guān)參數(shù),對計算模型進行敏感性分析,為跨臨界泄壓瞬態(tài)過程的準確計算提供參考。計算結(jié)果表明,加熱棒壁面是否發(fā)生溫度飛升取決于所選用的臨界熱流密度和最小膜態(tài)沸騰溫度的值;驟冷前沿模型的使用可以實現(xiàn)壁面再濕潤,降低壁面溫度。

跨臨界泄壓;瞬態(tài)傳熱;ATHLET 3.0程序;SWAMUP

超臨界水冷堆(Supercritical Water Cooled Reactor,SCWR)是六種第四代核反應(yīng)堆之一。相比常規(guī)壓水堆,超臨界水冷堆熱效率高,系統(tǒng)結(jié)構(gòu)簡化,具有更高的經(jīng)濟性和安全性[1-3]。其正常運行工況下系統(tǒng)壓力約為25MPa;在停堆或失水等事故工況下將經(jīng)歷跨臨界泄壓過程(即系統(tǒng)壓力由超臨界狀態(tài)降低到擬臨界點22.064MPa以下),而在臨界點附近,臨界熱流密度(CHF)的值很小,極易發(fā)生沸騰臨界,導(dǎo)致加熱壁面溫度急劇升高[4,5],因此該瞬態(tài)過程是超臨界水冷堆失水事故安全分析的關(guān)鍵。目前,跨臨界泄壓瞬態(tài)過程可以通過系統(tǒng)程序進行計算,但缺乏有效的實驗驗證[6]。

本文依托上海交通大學的超臨界流體多功能實驗回路(Supercritical WAter MUlti Purpose loop,SWAMUP)跨臨界泄壓過程的實驗[7],利用德國核安全中心(GRS)開發(fā)的系統(tǒng)程序ATHLET3.0進行建模計算[8],分析跨臨界泄壓過程傳熱特性。通過對次臨界區(qū)臨界熱流密度、最小膜態(tài)沸騰溫度、驟冷前沿模型等計算模型進行敏感性分析,研究跨臨界泄壓過程的傳熱原理,為實現(xiàn)該過程的準確計算提供方向。

1 SWAMUP實驗裝置介紹及ATHLET建模

作為歐盟第七框架協(xié)議資助項目“超臨界水冷堆燃料性能驗證實驗(SCWR-FQT)”的參與單位,上海交通大學的超臨界流體多功能實驗回路(SWAMUP)承擔對SCWR-FQT回路的安全許可申請的驗證工作之一,即跨臨界泄壓過程的實驗研究。SWAMUP實驗段結(jié)構(gòu)如圖1所示。2×2的直流加熱棒束通過繞絲定位,安裝在矩形流道中。經(jīng)過預(yù)熱后的冷卻水通過實驗段頂部進入第一流程,沿第一流程向下流動,同時通過矩形流道與第二流程進行熱量交換,吸收熱量。冷卻水在第一流程底部進行攪混后向上流動進入第二流道,冷卻加熱棒。實驗結(jié)果顯示,如圖2所示,當系統(tǒng)壓力跨臨界壓力時,靠近出口方向壁面溫度飛升,之后隨著壓力的進一步降低,壁面溫度回落。

圖1 實驗段示意圖Fig.1 Sketch of the test section

圖2 不同高度處的SWAMUP實驗結(jié)果Fig.2 SWAMUP measurement consequence at various axial positions

圖3是SWAMUP實驗段的ATHLET 3.0模型。將實驗段簡化為兩個平行管道,右側(cè)通道代表第一流程,左側(cè)通道代表第二流程,兩流程之間通過導(dǎo)熱構(gòu)件進行傳熱。加熱段長度為750mm,計算中將其劃分為100個控制體。根據(jù)SWAMUP跨臨界泄壓實驗,選取如表1所示瞬態(tài)工況進行計算。計算時,給定溫度、流量為入口邊界條件,出口邊界條件為壓力。超臨界泄壓瞬態(tài)過程中,入口溫度、入口流量、加熱棒功率、降壓速率等均保持恒定,系統(tǒng)壓力經(jīng)大約280s的時間從25MPa下降到17MPa。

圖3 實驗段ATHLET模型Fig.3 ATHLET model of test section

表1 跨臨界泄壓瞬態(tài)實驗相關(guān)參數(shù)Table 1 test parameters of pressure transient test case

2 瞬態(tài)計算結(jié)果

次臨界壓力下,當系統(tǒng)壓力接近擬臨界點22.064MPa時,臨界熱流密度(CHF)很低,即使在較低的熱流密度和高過冷度下,也很容易發(fā)生沸騰臨界,出現(xiàn)加熱棒壁面溫度急劇升高現(xiàn)象。對于跨臨界泄壓過程,為了確保安全分析計算的可靠性,理解該過程的傳熱機理顯得尤為重要。根據(jù)圖2所示實驗結(jié)果,該泄壓過程,主要關(guān)心以下三個方面的問題:

(1) 跨臨界壓力時,是否發(fā)生沸騰臨界,導(dǎo)致加熱棒壁溫飛升;

(2) 壁面所能達到的最高溫度;

(3) 壁溫升高后,隨著冷卻劑的不斷流入,壁面是否能夠得到再濕潤。

如圖4所示,是ATHLET3.0的計算結(jié)果。可見ATHLET3.0可以實現(xiàn)該跨臨界泄壓過程的計算,但計算結(jié)果與實驗結(jié)果(見圖2)差別很大。計算結(jié)果顯示,次臨界壓力下,未發(fā)生沸騰臨界,傳熱模式處在核態(tài)沸騰區(qū),傳熱系數(shù)較大,不會發(fā)生加熱棒壁面溫度飛升。該過程中,核態(tài)沸騰傳熱系數(shù)隨壓力降低而升高,因此加熱棒壁面溫度隨系統(tǒng)壓力降低而降低。

圖4 ATHLET3.0跨臨界泄壓瞬態(tài)計算結(jié)果Fig.4 the calculation result of the trans-critical transient in ATHLET 3.0

ATHLET3.0程序依據(jù)圖5所示沸騰曲線來選擇傳熱模式。先通過壁面溫度判斷換熱區(qū)間(單相區(qū)、核態(tài)沸騰區(qū)、過渡沸騰區(qū)、膜態(tài)沸騰區(qū)),進而通過平衡含氣率和空泡份額來選擇傳熱系數(shù)的計算公式。其中,臨界熱流密度溫度(TCHF)、最小膜態(tài)沸騰溫度(TMFB)、再濕潤溫度(TREW)、重返核態(tài)沸騰溫度(TRNB)對于傳熱模式選擇起著至關(guān)重要的作用[8]。在計算過程中,若以上判據(jù)溫度的選擇和傳熱系數(shù)的計算公式的選擇合理,便可以得到較準確的計算結(jié)果。本文選取其中的一些影響因素進行討論,分析超臨界泄壓瞬態(tài)計算對傳熱模型的敏感性。

圖5 ATHLET傳熱分區(qū)[9]Fig.5 ATHLET main heat transfer level

3 模型敏感性分析

3.1 臨界熱流密度的影響

ATHLET中,將臨界熱流密度通過公式轉(zhuǎn)換為臨界熱流密度溫度(TCHF)使用,作為是否發(fā)生沸騰臨界的判據(jù),若壁面溫度高于臨界熱流密度溫度(TCHF),傳熱模式便會從核態(tài)沸騰區(qū)跳轉(zhuǎn)到過渡沸騰區(qū)。類似地,由臨界熱流密度公式推導(dǎo)出重返核態(tài)沸騰溫度(TRNB),作為能否重新返回核態(tài)沸騰的判據(jù)。當沸騰臨界已經(jīng)發(fā)生,若壁面溫度低于TRNB,則傳熱模式將從過渡沸騰跳轉(zhuǎn)到核態(tài)沸騰區(qū)。然而,ATHLET3.0所提供的臨界熱流密度的計算公式,在高壓條件下預(yù)測CHF的準確性均較差。但計算時,用戶可以通過調(diào)節(jié)臨界熱流密度的系數(shù),同時實現(xiàn)臨界熱流密度溫度(TCHF)和重返核態(tài)沸騰溫度(TRNB)的調(diào)節(jié)。故本文通過調(diào)節(jié)臨界熱流密度系數(shù),改變計算過程中臨界熱流密度溫度和重返核態(tài)沸騰溫度的大小,觀察計算結(jié)果的變化。如圖6所示,選出三組典型的計算結(jié)果進行討論。計算中,所選取的臨界熱流密度計算公式為:ATHLET所提供的所有CHF計算公式的結(jié)果中取最小值。

圖6(a)所示結(jié)果對應(yīng)臨界熱流密度系數(shù)為0.2(即臨界熱流密度溫度和重返核態(tài)沸騰溫度的計算值均乘上0.2,下文提及最小膜態(tài)沸騰溫度等均與此類似)。相比于圖4所示結(jié)果,臨界熱流密度減小之后,靠近第二流程出口端壁面溫度在360s左右開始出現(xiàn)飛升。根據(jù)ATHLET計算結(jié)果輸出文件,隨系統(tǒng)壓力降低,臨界熱流密度溫度和重返核態(tài)沸騰溫度的計算值減小,因此靠近出口端加熱壁面溫度高于臨界熱流密度溫度(TCHF),發(fā)生沸騰臨界,傳熱模式由核態(tài)沸騰變?yōu)檫^渡沸騰進而變?yōu)槟B(tài)沸騰,傳熱系數(shù)減小。最終加熱壁面最高溫度為639.532℃。而靠近入口端壁面溫度尚未高于臨界熱流密度溫度,不發(fā)生沸騰臨界,因此并沒有發(fā)生壁面溫度飛升。

圖6(b)給出臨界熱流密度系數(shù)為0.1,即相比圖6(a)臨界熱流密度的值進一步減小。可以看出,第二流程所有的控制體加熱壁面溫度均高于臨界熱流密度溫度,發(fā)生沸騰臨界。傳熱模式經(jīng)過渡沸騰跳轉(zhuǎn)為膜態(tài)沸騰,傳熱系數(shù)減小,因此加熱棒壁面溫度急劇升高,降壓結(jié)束后,壁面溫度趨于穩(wěn)定值,最高溫度為639.532℃。

圖6 不同臨界熱流密度系數(shù)下的計算結(jié)果Fig.6 calculation results with different critical heat flux factor(a) CHF0.2; (b) CHF0.1; (c) CHF0.01

圖6(c)所示臨界熱流密度系數(shù)為0.01,第二流程最靠近出口的控制體在360s左右開始出現(xiàn)沸騰臨界。根據(jù)ATHLET計算結(jié)果輸出文件,隨系統(tǒng)壓力降低,臨界熱流密度溫度和重返核態(tài)沸騰溫度的計算值減小。起初壁面溫度高于臨界熱流密度溫度之后,傳熱模式變?yōu)檫^渡沸騰,之后因為壁面溫度繼續(xù)降低至低于重返核態(tài)沸騰溫度,傳熱模式再次變?yōu)楹藨B(tài)沸騰,如此傳熱模式在核態(tài)沸騰與過渡沸騰之間跳變,此時加熱棒壁面溫度并沒有升高,傳熱系數(shù)與圖3相比差別不大。直到400s左右,壁面溫度高于最小膜態(tài)沸騰溫度,傳熱模式轉(zhuǎn)變?yōu)槟B(tài)沸騰,傳熱系數(shù)迅速下降,壁面溫度開始急劇升高并趨于穩(wěn)定值,最高溫度為639.537℃。靠近入口方向上的控制體,沒有發(fā)生加熱棒壁面溫度飛升現(xiàn)象,傳熱模式一直在核態(tài)沸騰與過渡沸騰之間跳動。

以上三組計算結(jié)果可見,臨界熱流密度系數(shù)的大小影響傳熱模式,從而影響加熱棒壁面溫度。對于發(fā)生溫度飛升的壁面,其溫度將會趨于一個穩(wěn)定值而不再下降,且所能達到的最高溫度幾乎相同。

3.2 最小膜態(tài)沸騰溫度的影響

最小膜態(tài)沸騰溫度(TMFB)是過渡沸騰與膜態(tài)沸騰的分界,壁面溫度大于最小膜態(tài)沸騰溫度即進入膜態(tài)沸騰區(qū)。本節(jié)討論改變最小膜態(tài)沸騰溫度系數(shù)的計算結(jié)果(見圖7)。

圖7(a)所示為最小膜態(tài)沸騰溫度系數(shù)為0.9時的計算結(jié)果。400s左右,靠近出口端壁面溫度高于最小膜態(tài)沸騰溫度,傳熱模式由核態(tài)沸騰區(qū)直接轉(zhuǎn)變到膜態(tài)沸騰區(qū),傳熱系數(shù)減小,壁面溫度開始迅速升高。靠近入口端方向,因壁面溫度相對低,未超過最小膜態(tài)沸騰溫度,壁面溫度不發(fā)生突變。

圖7(b)最小膜態(tài)沸騰溫度系數(shù)為0.8。此時計算過程中所使用的最小膜態(tài)沸騰溫度更小,第二流程所有的控制體傳熱模式均將達到膜態(tài)沸騰,同時相繼出現(xiàn)壁面溫度飛升。

圖7 不同最小膜態(tài)沸騰溫度系數(shù)下的計算結(jié)果Fig.7 calculation results with different minimum film boiling temperature factor(a) TMFB0.9;(b) TMFB0.8

相比改變臨界熱流密度大小的計算結(jié)果,調(diào)節(jié)最小膜態(tài)沸騰溫度系數(shù),沸騰臨界后的換熱區(qū)直接從核態(tài)沸騰跳轉(zhuǎn)到膜態(tài)沸騰而不經(jīng)過渡沸騰區(qū)的過渡。計算結(jié)果對最小膜態(tài)沸騰溫度大小更敏感,系數(shù)調(diào)至0.9即帶來很大的不同。但兩種情況下得到的最高壁面溫度均為639.5℃,發(fā)生溫度飛升之后壁面不出現(xiàn)再濕潤現(xiàn)象。

3.3 驟冷前沿模型的影響

在ATHLET中,驟冷前沿(Quench front)模型定義了一種軸向的傳熱機制,使傳熱模式從膜態(tài)沸騰返回過渡沸騰或者核態(tài)沸騰,從而實現(xiàn)加熱壁面的再濕潤。在這種模型中,定義加熱壁面溫度高于Leidenfrost溫度的部分為干壁面,傳熱模式為膜態(tài)沸騰;類似地,定義加熱壁面溫度低于Leidenfrost溫度的部分為濕壁面,傳熱模式為過渡沸騰或者核態(tài)沸騰。干濕壁面的分界點稱為“驟冷前沿”。計算過程中,驟冷前沿的推進速度由Yamanouchi公式計算得到驟冷前沿模型作用下,會將驟冷前沿所到過的壁面溫度降到Leidenfrost溫度以下。[9]根據(jù)Hein的理論[10],驟冷前沿模型可以用于跨臨界泄壓過程的計算。據(jù)4.1節(jié)和4.2節(jié)的計算結(jié)果,臨界熱流密度系數(shù)為0.1和最小膜態(tài)沸騰溫度系數(shù)為0.8兩種條件下,均發(fā)生壁面溫度飛升,且沒有再濕潤現(xiàn)象發(fā)生,故在這兩種條件下,在壁面溫度達到極限值之后,添加驟冷前沿模型,所得計算結(jié)果如圖8所示。

圖8所示在510s添加驟冷前沿模型之后,傳熱模式由膜態(tài)沸騰變?yōu)楹藨B(tài)沸騰,傳熱系數(shù)增大,加熱棒壁面得到冷卻,從入口到出口壁面溫度依次回落。可見選擇在510s添加驟冷前沿模型,不影響壁面溫度的最高值,可以實現(xiàn)壁面再濕潤過程的模擬。但目前,尚且無法準確判斷添加驟冷前沿模型的時間。

圖8 驟冷前沿模型作用下的計算結(jié)果Fig.8 caculation with quench front model(a) CHF0.1; (b) TMFB0.8

4 結(jié)論

本文使用ATHLET3.0程序,對SWAMUP回路進行建模,實現(xiàn)跨臨界泄壓瞬態(tài)過程的計算,目前計算結(jié)果與實驗結(jié)果差別仍然較大,因此本文進行進一步計算得到不同臨界熱流密度系數(shù)及不同最小膜態(tài)沸騰溫度系數(shù)下,壁面溫度、傳熱模式、傳熱系數(shù)等的變化。可以得出以下結(jié)論:

(1) 跨臨界泄壓瞬態(tài)計算,受臨界熱流密度、最小膜態(tài)沸騰溫度的影響很大,這兩個參數(shù)的選擇將決定加熱棒壁面溫度是否會發(fā)生飛升;

(2) 沸騰臨界后的傳熱方式是過渡沸騰還是膜態(tài)沸騰影響壁面溫度計算結(jié)果,膜態(tài)沸騰傳熱系數(shù)更低,壁面溫度迅速升高;

(3) 驟冷前沿模型的加入,可以實現(xiàn)壁面再濕潤過程的計算,計算的可靠性還需要進一步驗證;

(4) 將來,需要對次臨界區(qū)臨界熱流密度和最小膜態(tài)沸騰溫度的值進行完善,對沸騰曲線的走向以及對再濕潤的條件做更多的研究。

致謝

感謝德國GRS核安全中心提供ATHLET程序。

[1] Oka Y. Review of High Temperature Water and Steam Cooled Reactor Concepts[J]. 2002.

[2] 程旭, 劉曉晶. 超臨界水冷堆國內(nèi)外研發(fā)現(xiàn)狀與趨勢[J]. 原子能科學技術(shù), 2008, 42(2):167-172.

[3] Oka Y, Koshizuka S, Ishiwatari Y, et al. Super Light Water Reactors and Super Fast Reactors: Supercritical-Pressure Light Water Cooled Reactors[M]. Springer Science & Business Media, 2010.

[4] Se-Young CHUN, Sung-Deok HONG, Hiroshige KIKURA, et al. Critical Heat Flux in a Heater Rod Bundle Cooled by R-134a Fluid near the Critical Pressure[J]. Journal of Nuclear Science and Technology, 2007, 44(9):1189-1198.

[5] Mawatari T, Mori H. An experimental study on characteristics of post-CHF heat transfer in the high subcritical pressure region near to the critical pressure[J]. Journal of Thermal Science & Technology, 2016, 11(1).

[6] Schulenberg T, Visser D C. Thermal-hydraulics and safety concepts of supercritical water cooled reactors[J]. Nuclear Engineering & Design, 2013, 264:231-237.

[7] 汪子迪, 曹臻, 劉曉晶,等. 超臨界水實驗回路熱工水力分析[J]. 原子能科學技術(shù), 2015(7):1191-1199.

[8] G. Lerchel, H.Austregesilo. ATHLET Mod3.0 Cycle A, User’s Manual, GRS-p-1/Vol.1,Rev.6.(2012).

[9] H.Austregesilo,C.Bals. ATHLET Mod3.0 Cycle A, Models and Methods, GRS-p-1/Vol. 4, Rev.3. (2014).

[10] K?hler W, Hein D. Influence of the wetting state of a heated surface on heat transfer and pressure loss in an evaporator tube[R]. Nuclear Regulatory Commission, Washington, DC (USA). Office of Nuclear Regulatory Research; Kraftwerk Union AG, Erlangen (Germany, FR), 1986.

SensitiveAnalysisofHeatTransferModelDuringTrans-criticalDepressurization

SONGMei-qi,LIUXiao-jing,CHENGXu

(School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240 , China)

During shutdown and loss of coolant accident conditions, Supercritical Water Cooled Reactor (SCWR) will undergo a pressure decrease from supercritical to subcritical pressure. And critical heat flux is likely to occur even at relatively low heat flux near critical pressure. Then, the heating wall will sharply heat up which causes rapid temperature increase. Hence, trans-critical transient are of crucial importance for LOCA safety analysis. Although the process can be calculated by several system codes, a reliable validation work is not available for now. In order to analyze the heat transfer mechanism, the trans-critical transient of Supercritical WAter MUltiPurpose loop (SWAMUP) with 2×2 rod bundle in Shanghai Jiao Tong University (SJTU) has been simulated by ATHLET 3.0, which is a kind of system code developed by GRS. And the calculations with different heat transfer parameters such as critical heat flux, minimum film boiling temperature and quench front model were carried out. The calculation results indicate that whether the heating wall temperature will sharply increase depends on the value of critical heat flux and minimum film boiling temperature. And the quench front model can bring rewetting process to the hot wall. In the future, the heat transfer model and the criterion for rewetting phenomena near critical pressure should be developed.

Trans-critical depressurization;Transient heat transfer;ATHLET 3.0;SWAMUP

2017-05-17

Science & Technology Commission of Shanghai Municipality (Grant No. 17580711400)

宋美琪 (1992—),女,山東滕州人,碩士研究生,現(xiàn)從事核科學與工程方面研究

TL331

A

0258-0918(2017)06-1053-08

主站蜘蛛池模板: 高清乱码精品福利在线视频| 精品久久香蕉国产线看观看gif | 国产美女免费| 国产成人福利在线| 欧美日韩国产成人高清视频| 在线国产欧美| 污网站免费在线观看| 国产无码性爱一区二区三区| 四虎国产在线观看| 亚洲一级色| 久久国产亚洲欧美日韩精品| 手机在线免费不卡一区二| 她的性爱视频| 国产一区二区三区免费观看| 中文天堂在线视频| 午夜限制老子影院888| 成人午夜福利视频| 97视频在线精品国自产拍| 亚洲一区二区三区国产精品| 99久久精品免费看国产电影| 国产黑丝一区| 黄色不卡视频| 免费在线a视频| 999国产精品| 久久香蕉欧美精品| 国产哺乳奶水91在线播放| 国产福利免费在线观看| 久久成人免费| 在线观看国产黄色| 在线另类稀缺国产呦| 国产成人AV综合久久| 久青草免费视频| 国产在线精品人成导航| 嫩草在线视频| 影音先锋亚洲无码| 亚洲视频免费在线看| 亚洲成人免费看| 亚洲国产91人成在线| 欧洲熟妇精品视频| 国产在线观看精品| 91久久精品日日躁夜夜躁欧美| 国产成人区在线观看视频| 欧美日韩国产精品va| 五月婷婷导航| 婷婷亚洲最大| 日韩欧美中文字幕在线韩免费| 亚洲视频三级| 国产成人亚洲毛片| 天天综合色网| 99re在线视频观看| 精品无码视频在线观看| 美女免费精品高清毛片在线视| 久久99精品久久久久纯品| 欧美激情成人网| 波多野结衣视频一区二区| 五月天久久婷婷| 久久久久久久蜜桃| 日韩无码白| 国产激情在线视频| 国产成人亚洲精品无码电影| 看国产毛片| 日韩毛片基地| 欧美狠狠干| 亚洲二区视频| 国产交换配偶在线视频| 99久久精品国产精品亚洲| 9久久伊人精品综合| 久久国产成人精品国产成人亚洲 | 激情亚洲天堂| 国产精品香蕉| 国产精品丝袜视频| 91久久国产热精品免费| 亚洲色图欧美一区| 婷婷在线网站| 国产人成乱码视频免费观看| 久久久久久久久18禁秘 | 中文字幕人妻无码系列第三区| 精品人妻无码中字系列| 亚洲欧美成人网| 国产青青操| 国产大片喷水在线在线视频| 国产麻豆va精品视频|