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

大功率先進(jìn)壓水堆IVR有效性評價中熔池?fù)Q熱研究

2014-05-25 00:33:40劉曉晶楊燕華
原子能科學(xué)技術(shù) 2014年2期
關(guān)鍵詞:模型

鮑 晗,金 越,劉曉晶,楊燕華

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

大功率先進(jìn)壓水堆IVR有效性評價中熔池?fù)Q熱研究

鮑 晗,金 越,劉曉晶,楊燕華

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

熔融物堆內(nèi)滯留-壓力容器外部冷卻(IVR-ERVC)是一種重要的核電廠嚴(yán)重事故緩解措施。當(dāng)前針對IVR有效性評價的方法主要是基于集總參數(shù)模型對下封頭熔池?fù)Q熱進(jìn)行分析。在大功率先進(jìn)壓水堆熔池集總參數(shù)法計算中,堆芯重量變大、壓力容器尺寸增加會導(dǎo)致熔池自然對流換熱中的瑞利數(shù)Ra′增大。通過使用集總參數(shù)分析程序,對比研究熔池氧化層各換熱模型的適用范圍,計算大功率先進(jìn)壓水堆高瑞利數(shù)條件下穩(wěn)態(tài)熔池的自然對流換熱,模擬兩層穩(wěn)態(tài)熔池模型中壓力容器外壁面的熱流密度分布,對其進(jìn)行選定嚴(yán)重事故序列下的IVR-ERVC有效性評價,并對堆內(nèi)構(gòu)件設(shè)計提出建議。

熔融物堆內(nèi)滯留-壓力容器外部冷卻;大功率先進(jìn)壓水堆;集總參數(shù)模型;自然對流換熱

熔融物堆內(nèi)滯留-壓力容器外部冷卻(IVRERVC)是一種重要的核電廠嚴(yán)重事故緩解措施。IVR熱負(fù)荷失效是指由于熔池內(nèi)衰變熱引起壓力容器壁面的熱流密度超過其外部自然循環(huán)的臨界熱流密度(CHF),造成壁面熔穿,放射性物質(zhì)泄漏至壓力容器外部。當(dāng)下封頭熔池達(dá)到最終穩(wěn)定熱力狀態(tài)時,崩塌至下封頭的堆芯燃料碎片體積達(dá)到最大,熔池溫度達(dá)到最高,金屬層上表面向上熱輻射路徑的熱阻達(dá)到最大,此時壓力容器下封頭壁面的熱負(fù)荷達(dá)到最大值[1]。因此,對于下封頭穩(wěn)態(tài)熔池?zé)崃μ匦缘难芯渴窃u價反應(yīng)堆IVR熱負(fù)荷失效的重點。

目前,在役的核電站和先進(jìn)堆型大都采用IVR緩解措施并對其有效性進(jìn)行了評價。在AP600、AP1000的穩(wěn)態(tài)熔池?zé)崃μ匦匝芯恐校饕獞?yīng)用的均是集總參數(shù)法[1-2]。目前,國際上正針對更大功率的核反應(yīng)堆進(jìn)行IVR-ERVC的論證和設(shè)計工作,如韓國的APR1400[3]。在針對大功率先進(jìn)壓水堆(1 400/1 700MW級)的熔池計算中,由于堆芯重量變大,熱功率提升,壓力容器尺寸增加,導(dǎo)致瑞利數(shù)Ra′增大,因此需要選取適合高瑞利數(shù)的換熱關(guān)系式進(jìn)行模擬計算,這對進(jìn)行大功率先進(jìn)壓水堆IVRERVC有效性研究有著重要意義。

本文通過分析常用氧化池自然對流換熱關(guān)系式及其適用范圍,并使用集總參數(shù)分析程序模擬在選定嚴(yán)重事故下的熔池自然對流換熱,以研究不同瑞利數(shù)條件下氧化池自然對流換熱關(guān)系式對大功率先進(jìn)壓水堆熔池?zé)崃μ匦缘挠绊懀瑫r對其IVR-ERVC的有效性進(jìn)行評價,并對堆內(nèi)構(gòu)件設(shè)計提出建議。

1 集總參數(shù)分析模型

本研究所開發(fā)的程序采用集總參數(shù)分析

模型,根據(jù)熔池內(nèi)熱平衡的熔池分層傳熱模型計算壓力容器下封頭外壁面局部熱流密度,與局部CHF比較后得到IVR熱負(fù)荷失效的評價結(jié)果。

1.1 穩(wěn)態(tài)熔池分層模型

受實驗條件和經(jīng)驗的限制,現(xiàn)在還沒有一確定的實驗或方法可確定發(fā)生堆芯熔化的嚴(yán)重事故后下封頭內(nèi)的熔池結(jié)構(gòu)。美國加州大學(xué)圣巴巴拉分校(UCSB)對AP600進(jìn)行的熔池傳熱計算時采用了兩層模型[1]。這種結(jié)構(gòu)包絡(luò)整個中間狀態(tài)的熱負(fù)荷狀況,并代表任何IVR事故序列均會出現(xiàn)的最終狀態(tài),包括從初始熔融到所有堆芯碎片均重定位至下封頭的最終穩(wěn)態(tài),具有一定的保守性。之后美國核管會(NRC)評價AP1000和INEEL評價韓國APR1400均提及這種保守包絡(luò)結(jié)構(gòu)。

由于UCSB對AP600進(jìn)行的IVR評價已得到NRC的承認(rèn)并進(jìn)行了完整的驗證[4],因此本研究中所開發(fā)的程序最后確定采用UCSB的兩層模型,將堆芯熔化后壓力容器下封頭內(nèi)的穩(wěn)態(tài)熔池模型由上至下分為金屬層和氧化池,其中氧化池邊界由于遇冷凝結(jié)形成一硬殼,熔池模型示意圖示于圖1。下層氧化池是熱源區(qū),內(nèi)部進(jìn)行自然對流換熱。氧化池中的堆芯熔融物成分主要由熔化的堆芯燃料UO2和燃料包殼中的鋯被氧化后形成的ZrO2構(gòu)成;金屬層不含內(nèi)熱源,主要由熔化的燃料棒包殼、未氧化的Zr合金、堆內(nèi)構(gòu)件中熔化的不銹鋼及其他堆芯結(jié)構(gòu)熔化材料等組成。

圖1 穩(wěn)態(tài)熔池兩層模型示意圖Fig.1 Schematic of two-layered stratified molten pool configuration

1.2 熱平衡方程

兩層模型采用的熱平衡方程如下。

金屬層:

氧化池:

其中:Q?為單位衰變熱功率,W/m3;q″為熱流密度,W/m2;V為體積,m3;A為換熱面積,m2;下標(biāo)l表示金屬層,o表示氧化池,w表示壓力容器壁面,cr表示氧化硬殼,b表示底部,t表示頂部,s表示側(cè)面,i表示內(nèi)壁面。

由于氧化池上層硬殼及金屬層很薄,在氧化池上部硬殼的體積忽略不計。

1.3 熔池?fù)Q熱分層計算

1)氧化池

根據(jù)選定模型的換熱關(guān)系式和衰變熱數(shù)據(jù)可得到氧化池向上平均熱流密度和向下平均熱流密度,進(jìn)而計算溫度、硬殼厚度、局部壓力容器壁厚和熱流密度等參數(shù)。

氧化池向上與向下平均努賽爾數(shù)Nu比為:

則氧化池向下、向上平均熱流密度為:

氧化池向下局部熱流密度為:

硬殼厚度δcr,s(θ)的計算公式為:

壓力容器內(nèi)壁面熱流密度為:

壓力容器壁面厚度的計算公式為:

同時有壓力容器外壁面熱流密度為:

其中:θ為角度;To,m為氧化硬殼熔化溫度,K;κcr為硬殼熱導(dǎo)率,W/(m·K);κw為壓力容器壁面熱導(dǎo)率,W/(m·K);Cboil為堆腔注水的核態(tài)沸騰系數(shù);Tsat為堆腔注水的飽和溫度,K。

2)金屬層

由熱平衡方程可得金屬層底部熱流密度:

其中:hl,b為金屬層下表面與氧化池上部硬殼間的換熱系數(shù),W/(m2·K);Tl,bulk為金屬層中心溫度,K;Tl,b為金屬層下表面溫度,K。

金屬層對其上表面的換熱關(guān)系式為:

其中:hl,t為金屬層上表面的換熱系數(shù),W/(m2· K);Tl,t為金屬層上表面溫度,K。

金屬層上表面對壓力容器上腔室大空間的輻射換熱計算公式為:

其中:Ts為上腔室構(gòu)件平均溫度,K;As為上腔室構(gòu)件表面積,m2;σ為Stefan-Boltzman常數(shù);εt為金屬層上表面發(fā)射率;εs為上腔室構(gòu)件發(fā)射率。

考慮到上腔室構(gòu)件的內(nèi)、外溫度不同,通過迭代計算出溫度。

金屬層對側(cè)面壓力容器壁面的換熱為:

其中,hl,s為金屬層側(cè)面與壓力容器壁面間的換熱系數(shù),W/(m2·K)。

壓力容器壁面厚度δw(θ)和壓力容器外壁面至外部流道的熱流密度可分別用式(11)、(12)求出。

3)金屬層上部空間輻射換熱模型

對金屬層上表面與上腔室構(gòu)件的換熱,模型考慮了金屬層上表面與上腔室構(gòu)件之間的熱輻射,上腔室構(gòu)件內(nèi)、外壁面之間的熱傳導(dǎo),以及上腔室構(gòu)件外壁面與壓力容器內(nèi)壁面之間的熱輻射,假設(shè)了初始溫度并通過迭代計算得到金屬層各表面和上腔室構(gòu)件的內(nèi)、外壁面溫度。

在式(15)中,假設(shè)上腔室構(gòu)件的平均溫度為Ts,若考慮到上腔室構(gòu)件的內(nèi)、外溫度不同,則需通過迭代計算出溫度。用Ts,i、Ts,o分別表示上腔室構(gòu)件內(nèi)、外壁面溫度,則有:

1.4 衰變熱模型

在兩層模型中,假設(shè)只有下層氧化池內(nèi)有衰變熱,上層金屬層內(nèi)部無衰變熱。

式中,Pdecay-tot為熔池總衰變熱,W。

1.5 臨界熱流密度模型

實際RPV外壁面CHF數(shù)據(jù)應(yīng)根據(jù)針對性的CHF實驗獲得。由于目前大功率先進(jìn)壓水堆的CHF實驗尚未完成,本工作用于計算1 700MW級的CHF關(guān)系式根據(jù)AP600的關(guān)系式類推得到:

2 氧化池自然對流換熱模型

在堆熔池分析中,氧化池的流動屬于含內(nèi)熱源的半球形區(qū)域自然對流,金屬層的流動屬于不含內(nèi)熱源的圓形平板間自然對流。

在針對大功率先進(jìn)壓水堆(如1 700MW級)熔池集總參數(shù)法計算中,由于堆芯重量變大,熱功率提升,壓力容器尺寸增加,導(dǎo)致Ra′增大,需選取適合高瑞利數(shù)的換熱關(guān)系式進(jìn)行模擬。Ra′與熔池氧化層衰變功率Q和氧化池高度H成正比,基于Theofanous分析AP600的IVR有效性驗證計算時所得到的熔池氧化層Ra′分布[1],由熱功率和熔池尺寸比例簡單類推出1 700MW級反應(yīng)堆熔池氧化層Ra′約為1017量級。表1列出常用的氧化池向上、向下對流換熱關(guān)系式及其適用范圍[2]。

表1 氧化池自然對流換熱關(guān)系式及適用范圍Table 1 Heat transfer correlations for ceramic pool

以上兩種關(guān)系式可計算出氧化池衰變熱向上和向下的分配比例,向下傳遞的衰變熱在傳給壓力容器壁面時是隨角度分布的。UCLA實驗提出了局部換熱系數(shù)與整體換熱系數(shù)比值隨角度分布的關(guān)系式[5]:

其中,φ為氧化池液面所達(dá)到的角度,°。

mini-ACOPO實驗也總結(jié)出類似的熱流密度隨角度分布的曲線[6]:式中,θp為氧化池液面所達(dá)到的角度。

圖2示出氧化池向上、下對流換熱關(guān)系式與ACOPO實驗數(shù)據(jù)的對比,在超過關(guān)系式適用范圍時使用外推方法進(jìn)行模擬。圖3示出氧化池向下半球局部熱流密度隨角度的變化。通過對比可發(fā)現(xiàn),ACOPO實驗得出的關(guān)系式在氧化池的Ra′范圍內(nèi)與實驗數(shù)據(jù)的偏差相對較小,且適用范圍更廣,Ra′可達(dá)2× 1016,最接近1 700MW級等大功率先進(jìn)壓水堆簡單推導(dǎo)得到的范圍。

圖2 氧化池向上(a)、向下(b)對流換熱關(guān)系式與實驗數(shù)據(jù)的對比Fig.2 ACOPO data vs.upward(a)and downward(b)heat transfer correlations

圖3 氧化池向下半球局部熱流密度隨角度的變化Fig.3 Downward heat flux vs.angular

3 計算結(jié)果討論

使用熔池集總參數(shù)分析程序,從氧化池至金屬層再到上腔室空間逐層計算各層傳熱。同時假設(shè)在堆內(nèi)熔融物坍塌至壓力容器下封頭之前堆腔中已注滿水,之后堆腔注水系統(tǒng)的補水量大于堆腔中水的蒸發(fā)量,溫度恒定。最后得到沿壓力容器下封頭軸向各角度處的熱流密度分布并與臨界熱流密度比較,判斷壓力容器發(fā)生熱負(fù)荷失效的可能性。重點分析不同氧化池自然對流換熱模型對計算結(jié)果的影響。

3.1 事故序列選取

由于對大功率先進(jìn)壓水堆設(shè)計和嚴(yán)重事故序列不確定性的研究還處于初步階段,本文僅選取具有包絡(luò)性的大破口失水事故(熱管段雙端斷裂)作為評價大功率先進(jìn)壓水堆IVR有效性的事故序列。由于目前國內(nèi)外1 700MW級電站的預(yù)研剛剛開始,部分設(shè)計參數(shù)還未確定,無法使用MAAP等嚴(yán)重事故分析程序進(jìn)行熔池計算。考慮到1 700MW級電站和AP1000均采用非能動技術(shù),同時假設(shè)兩者有相似的堆芯燃料組件結(jié)構(gòu),用MAAP程序計算得到該事故發(fā)生后1 000MW級堆內(nèi)熔池中各不確定輸入?yún)?shù)后,通過熱功率比類推得到相同事故時間下1 700MW級堆內(nèi)熔池的相應(yīng)參數(shù),結(jié)果列于表2。

表2 堆內(nèi)熔池各不確定輸入?yún)?shù)Table 2 Uncertainties for selected severe accident

3.2 不同氧化池自然對流換熱模型的影響

將氧化池自然對流換熱關(guān)系式分為3組,分別為[2]:1)K-M模型(分別使用Kulacki-Emara的向上換熱關(guān)系式和Mayinger的向下?lián)Q熱關(guān)系式);2)mini-ACOPO模型(使用mini-ACOPO向上和向下?lián)Q熱關(guān)系式);3)ACOPO模型(使用ACOPO向上和向下?lián)Q熱關(guān)系式)。

根據(jù)3.1節(jié)選定的事故序列,模擬1 700MW級堆穩(wěn)態(tài)熔池壓力容器外壁面熱流密度,結(jié)果示于圖4。

圖4 不同換熱模型的壓力容器外壁面熱流密度Fig.4 RPV wall heat flux with different heat transfer models

由于沒有相關(guān)實驗或模擬計算作為驗證對比,本文僅以模型的保守性作為模型對熔池自然對流換熱和IVR-ERVC有效性影響的衡量標(biāo)準(zhǔn)。由圖4可看出,在兩層熔池模型軸向角度0°~39°內(nèi)(氧化池的下部區(qū)域),K-M模型最保守;在40°~73°(氧化池上部至金屬層底部區(qū)域),mini-ACOPO模型最為保守;在74°~86°(金屬層區(qū)域),ACOPO模型最為保守。考慮到在氧化池區(qū)域,壓力容器外壁面熱流密度與CHF相差較大,基本不會發(fā)生熱負(fù)荷失效,而在金屬層區(qū)域,熱流密度與CHF較為接近,甚至在氧化池與金屬層交界區(qū)域超過了局部CHF,因此在計算中可考慮使用在金屬層區(qū)域保守度最高的ACOPO模型進(jìn)行模擬計算。UCSB針對AP600的模型中使用的是mini-ACOPO的換熱關(guān)系式[1],本工作在計算1 700MW級反應(yīng)堆模型時將采用ACOPO換熱關(guān)系式進(jìn)行計算,仍使用UCSB的兩層模型結(jié)構(gòu)和分析方法。

同時考慮到ACOPO實驗所提出的經(jīng)驗關(guān)系式的瑞利數(shù)適用范圍與1 700MW級反應(yīng)堆推導(dǎo)得出的范圍更為接近,在1 700MW級堆熔池集總參數(shù)法中,可考慮選取ACOPO實驗得出的向上、向下?lián)Q熱關(guān)系式和mini-ACOPO的向下半球熱流密度分布曲線。

3.3 大功率先進(jìn)壓水堆IVR-ERVC有效性評價

不考慮應(yīng)力失效,當(dāng)壓力容器外壁面任一點的熱流密度超過臨界熱流密度時,即發(fā)生熱負(fù)荷失效。熱負(fù)荷失效是評價IVR有效性的唯一準(zhǔn)則。對具有包絡(luò)性的大破口失水事故(熱管段雙端斷裂),使用ACOPO實驗總結(jié)的氧化池自然對流換熱關(guān)系式,用集總參數(shù)分析程序計算得到的壓力容器外壁面熱流密度、壓力容器壁面厚度和氧化池側(cè)面硬殼厚度隨角度的變化,結(jié)果示于圖5~7。

圖5 壓力容器外壁面熱流密度隨角度的變化Fig.5 RPV wall heat flux vs.angular

圖6 壓力容器壁面厚度隨角度的變化Fig.6 RPV wall thickness vs.angular

圖7 氧化池側(cè)面硬殼厚度隨角度的變化Fig.7 Crust thickness vs.angular

由圖5可知,以類推得到的1 700MW級堆內(nèi)構(gòu)件設(shè)計,在熱段雙端斷裂的大破口失水事故中,堆芯熔化崩塌至下封頭形成穩(wěn)定熔池后,壓力容器外壁面熱流密度在氧化池與金屬層交界處(約73°)超過CHF,有發(fā)生熱負(fù)荷失效的可能。由圖6可知,壓力容器壁隨角度增加而逐步變薄,更易發(fā)生熔穿。在圖7中,金屬層側(cè)面無硬殼,故硬殼厚度為0。

3.4 金屬層熱聚集效應(yīng)的影響

由圖5可明顯看出,熱流密度在交界處發(fā)生顯著提升,導(dǎo)致這種情況出現(xiàn)的重要原因是金屬層的熱聚集效應(yīng)。當(dāng)堆芯熔化并塌落至下封頭形成分層熔池后,下封頭壁面受熱的熱流密度沿其軸向是非均勻分布的;上層金屬層由于質(zhì)量小、厚度薄,造成其熱流密度很大,形成熱聚集效應(yīng)。因此,建議在進(jìn)行反應(yīng)堆設(shè)計時需添加在堆內(nèi)構(gòu)件處的犧牲材料(如不銹鋼等),在堆芯熔化形成下封頭穩(wěn)定熔池時,會使上層金屬層厚度增加以緩解熱聚集效應(yīng)的影響。

將表2中1 700MW級堆熔池中的不銹鋼質(zhì)量稍微增加來模擬添加犧牲材料后落入下封頭熔池的結(jié)果,將程序輸入?yún)?shù)中的“熔池內(nèi)不銹鋼質(zhì)量”分別設(shè)為62、65和68t。圖8示出添加犧牲材料后壓力容器外壁面局部熱流密度與CHF的比值。由圖8可看出,隨不銹鋼質(zhì)量的增加,金屬層部分熱流密度比值越來越小,并低于1,即外壁面局部熱流密度小于CHF,不會發(fā)生熱負(fù)荷失效。

圖8 添加犧牲材料后壓力容器外壁面熱流密度與CHF的比值Fig.8 Ratio of RPV wall heat flux and CHF with adding sacrificial material

4 結(jié)論

本工作分析并提出了目前較適用于高瑞利數(shù)條件下熔池計算的自然對流換熱關(guān)系式,使用集總參數(shù)分析程序?qū)Υ蠊β氏冗M(jìn)壓水堆在高瑞利數(shù)條件下穩(wěn)態(tài)熔池的自然對流換熱進(jìn)行了計算,模擬了兩層穩(wěn)態(tài)熔池模型中壓力容器外壁面的熱流密度分布。結(jié)果表明:在包絡(luò)性的大破口失水事故(熱管段雙端斷裂)下,根據(jù)對大功率先進(jìn)壓水堆(1 700MW級)的初步設(shè)計,IVR-ERVC因金屬層熱聚集效應(yīng)有發(fā)生熱負(fù)荷失效的可能,但通過在堆內(nèi)構(gòu)件中添加犧牲材料(如不銹鋼等)會對金屬層的熱聚集效應(yīng)產(chǎn)生緩解作用,這對評價大功率先進(jìn)壓水堆在嚴(yán)重事故下的IVR-ERVC有效性具有重要意義。

[1] THEOFANOUS T G,LIU C,ADDITON S,et al.In-vessel coolability and retention of a core melt,DOE/ID-10460[R].USA:DOE,1996.

[2] Analysis of in-vessel retention and ex-vessel fuel coolant interaction for AP1000NUREG/CR-6849[R].USA:NRC,2004.

[3] REMPE J L,SUH K Y,CHEUNG F B,et al.In-vessel retention strategy for high power reactors,INEEL/EXT-04-02561[R].USA:INEEL,2005.

[4] REMPE J L,KNUDSON D L,ALLISON C M,et al.Potential for AP600in-vessel retention through ex-vessel flooding,INEEL/EXT-97-00779[R].USA:INEEL,1997.

[5] ASFIA F J,DHIR V K.An experimental study of natural convection in a volumetrically heated spherical pool bounded on top with a rigid wall[J].Nuclear Engineering and Design,1996,163:333-348.

[6] THEOFANOUS T G,LIU C,SCOTT J,et al.Natural convection in hemispherical enclosures at internal Rayleigh numbers up to 7×1014[R].Santa Barbara,USA:Center for Risk Studies and Safety University of California,1996.

Study of Natural Convection Heat Transfer in Molten Pool for Advanced Large Size PWR

BAO Han,JIN Yue,LIU Xiao-jing,YANG Yan-hua
(School of Nuclear Science and Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

In-vessel retention-external reactor vessel cooling(IVR-ERVC)is a key severe accident management strategy,and the present assessment method on the effectiveness of IVR is mainly based upon the lumped parameter model.With the enlargement of core mass and the RPV size,the increase of Rayleigh number was found in the calculation of nature convection heat transfer for advanced large size PWR.The applicability of several natural convection heat transfer correlations was studied and compared,and natural convection heat transfer under the high Rayleigh number was simulated by using lumped parameter analysis code.A two-layered stratified molten pool configuration was assumed to be developed in the reactor vessel lower head.The heat flow distribution at the outer boundary of RPV wall was simulated,and the code was applied to assess the feasibility of the IVR-ERVC strategy under a selected severe accident for advanced large size PWR.Recommendation for the designs of reactor core internals was proposed.

IVR-ERVC;advanced large size PWR;lumped parameter model;natural convection heat transfer

TL334

A

1000-6931(2014)02-0234-07

10.7538/yzk.2014.48.02.0234

2012-11-27;

2013-01-14

鮑 晗(1988—),男,湖北襄陽人,碩士研究生,核能科學(xué)與工程專業(yè)

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产一区二区精品福利 | 综合人妻久久一区二区精品| 色爽网免费视频| 欧美激情首页| 天天做天天爱夜夜爽毛片毛片| 日韩欧美国产另类| а∨天堂一区中文字幕| 国产一区二区免费播放| 免费jizz在线播放| 久久伊伊香蕉综合精品| 一级片免费网站| 欧美视频二区| 亚洲色婷婷一区二区| 亚洲一区波多野结衣二区三区| 日韩少妇激情一区二区| 日本91在线| 久久精品丝袜| 亚洲男女天堂| 69免费在线视频| 亚洲国产高清精品线久久| 日韩专区欧美| 日韩欧美91| 欧美啪啪视频免码| 无码丝袜人妻| 91亚洲精品国产自在现线| 精品国产女同疯狂摩擦2| 激情无码字幕综合| 亚洲精品动漫在线观看| 亚洲天堂视频在线观看| 不卡网亚洲无码| 九九九精品视频| 波多野结衣无码中文字幕在线观看一区二区 | 国产亚洲精品yxsp| 亚洲一级毛片在线观| 99视频精品全国免费品| 99成人在线观看| 91九色国产在线| 国产精品女主播| 十八禁美女裸体网站| 欧美人与牲动交a欧美精品| 国产97视频在线| 丰满人妻被猛烈进入无码| 欧美精品xx| 国产精品无码制服丝袜| 欧美有码在线| 亚洲丝袜第一页| 日韩无码黄色| 亚洲欧洲日韩国产综合在线二区| 成人小视频在线观看免费| 亚洲天堂视频在线观看| 亚洲综合日韩精品| 国产精品3p视频| 91偷拍一区| 日韩免费中文字幕| 国产视频自拍一区| 婷婷综合在线观看丁香| 色香蕉网站| 青草视频免费在线观看| 一级毛片在线播放免费观看| 日韩精品毛片人妻AV不卡| 狂欢视频在线观看不卡| 色综合久久无码网| 国产欧美在线观看一区| 日韩精品免费一线在线观看| 青青热久免费精品视频6| 1769国产精品视频免费观看| 亚洲美女一区二区三区| 3344在线观看无码| 波多野结衣久久精品| 国产成人区在线观看视频| 18黑白丝水手服自慰喷水网站| 国产精品思思热在线| 波多野吉衣一区二区三区av| 国产SUV精品一区二区6| 91视频首页| 999精品色在线观看| 国产成人综合在线视频| 99在线视频免费观看| 久久永久免费人妻精品| 日韩无码白| 亚洲无码高清视频在线观看| 国产日韩丝袜一二三区|