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

水-氫交換氫同位素體系HD/H2O、DT/D2O和HT/H2O分離性能模擬研究

2016-09-09 07:06:03胡石林尹玉國黃登高竇勤成
核化學與放射化學 2016年4期
關鍵詞:理論體系研究

吳 棟,阮 皓,胡石林,尹玉國,張 麗,黃登高,竇勤成

中國原子能科學研究院 特種材料專項工程部,北京 102413

?

水-氫交換氫同位素體系HD/H2O、DT/D2O和HT/H2O分離性能模擬研究

吳棟,阮皓,胡石林,尹玉國,張麗,黃登高,竇勤成

中國原子能科學研究院 特種材料專項工程部,北京102413

計算模擬應用于氫同位素分離領域,能夠方便、快捷地進行工藝條件分析。本工作采用數值模擬的方法對比研究了水-氫催化交換過程中HD/H2O、DT/D2O和HT/H2O三種氫同位素體系的分離性能。研究表明:在一定工藝條件下,三種體系均在操作溫度為343 K時達到最大的分離效果;隨著氣液比從1.0增大到3.0,最優操作溫度均從343 K降低到323 K,但是在此過程中,HT/H2O體系的分離效果受溫度的影響較小一些;在達到最大分離效果的目標下,HT/H2O體系需要的理論塔板數比HD/H2O和DT/D2O體系少,同時,在優化的工藝條件下,三體系氣相中氫同位素濃度在交換柱內分布曲線存在一定的差異。

水-氫交換;氫同位素分離;數值模擬;工藝條件

在氫同位素分離領域中,采用水-氫同位素催化交換的方法實現不同氫同位素之間分離的研究越來越廣泛。計算機模擬研究作為與試驗研究優勢互補的手段越來越受到科研工作者的重視,該方法在水-氫同位素催化交換工藝中的應用也越來越廣泛。近年來,基于平衡級模型的計算機模擬研究取得了一些成果。葉林森等[1]對填料塔內的H2/HDO體系進行模擬研究,研究了溫度、氣液比等對分離效果的影響;日本、俄羅斯和韓國等國家依據平衡級模型開發了相關的計算程序[2-4]。除此之外,夏修龍等[5]建立了動態模型,分析了 HD/H2O和HT/H2O兩體系催化交換工藝中電解池滯留量、進料位置和傳質系數對分離性能的影響。盡管如此,由于實驗條件的限制,尤其是T操作過程對輻射防護有較高的要求,至今未開展過水氫交換過程中氕(H)、氘(D)和氚(T)三種氫同位素的分離性能差別的實驗研究。因此,本工作依據平衡級模型,采用計算模擬的方法對HD/H2O、DT/D2O和HT/H2O三種體系氫同位素分離性能進行對比研究,研究目的為:(1) 系統探索溫度、氣液比和理論塔板數等工藝參數對三種體系分離效果的影響;(2) 針對不同分離體系獲得優化的工藝條件;(3) 分析三種體系分離性能的差異,獲得一般規律,以為工業應用提供有力的技術支撐。

1 模型與計算方法

1.1模型介紹

平衡級模型在水-氫交換過程中研究比較成熟,應用也比較廣泛,基本關系為物料平衡、反應平衡和氣液平衡關系。依據平衡關系建立數學方程,然后通過求解數學方程獲得描述催化交換過程的工藝參數。在液相催化交換工藝過程中,氫氣與液態水互相接觸實現不同氫同位素之間的分離,可以應用于H、D和T三種氫同位素之間的分離。該模型以含H和D兩種氫同位素的 HD/H2O體系為例進行說明[6]。在平衡級模型中,從交換柱底部進入的HD(H2)氣體和從頂部進入的液態H2O(HDO)發生氫同位素之間的交換反應,如圖1方程(a)所示。將交換柱內的兩步反應進行模型化處理:方程(b)表示的氣相反應過程發生在催化劑層;方程(c)表示的相交換過程發生在填料層。此時,氫氣通過填料層時氫氣中氘濃度不發生變化,液態水通過催化劑層時水中氘濃度也不發生改變。每一個催化劑層和填料層組合構成一個平衡級,也稱作理論塔板,每一個塔板上發生 (b)和(c)兩步反應,假設催化交換柱具有N級理論塔板。

圖1 水-氫交換催化交換模型圖Fig.1 Schematic model of liquid-phase catalytic exchange column

在第N級的催化劑層中,氘的物料平衡關系:

(1)

其中G和V分別是氫氣和水蒸氣的流量,y和z分別是氫氣和水蒸氣中的氘濃度,*代表平衡狀態。

其中在氘濃度比較低的情況下,氘在氫氣與水蒸氣的分離因子αg等于交換反應的平衡常數Kg,則有:

(2)

同樣,發生在填料床層的第N級氘的平衡關系:

(3)

其中L是液態水的流量,x是液態水中的氘濃度。

氘在液態水和水蒸氣的分離因子為:

(4)

聯立方程(1)—(4)可得:

(5)

其中:γg=G/V;γv=L/V

求解方程(5)可以得到第N級氫氣中氘濃度:

(6)

其中:

催化交換柱頂部和底部液態水中的氘濃度為:

(7)

從交換柱底部的質量守恒可以獲得底部水蒸氣中的D濃度:

(8)

在已知進料氫氣、液態水中的氘濃度和兩個分離因子的情況下,依據上述公式可以獲得交換柱頂部氫氣和底部液態水中的氘濃度。上述模型應用于DT/D2O和HT/H2O體系,總反應公式如方程(9)和(10)所示。

(9)

(10)

1.2分離因子

從計算模型中可以看出,分離因子的獲取在模型求解過程中至關重要。由于氕氘氚的同位素效應比較大,水-氫交換反應的平衡常數不等于1。通常采用分離因子來表征同位素效應的大小。對于水-氫同位素催化交換過程,在溫度比較高和氫同位素濃度比較低(一般小于8%)的情況下,分離因子和平衡常數看作一致。水-氫同位素的平衡常數由氣相反應和汽液相變兩步交換的平衡常數決定,并且前者只與溫度相關,而后者與溫度和氫同位素濃度兩者相關。

氣相反應的平衡常數一般采用實驗測定,然后通過對獲得的平衡常數值與溫度進行數據擬合,以獲得隨溫度變化的關系式并作為經驗公式使用。加拿大和俄羅斯的研究人員分別建立了針對不同氫同位素體系的經驗公式,這些經驗公式在實際中也有比較廣泛的應用[7]。除此之外,氣相反應的平衡常數還可以通過量子化學計算獲得。鐘正坤等[8]采用量子化學從頭計算方法,通過計算獲得各同位素分子的配分函數值,然后推導出HD/H2O和DT/D2O兩體系的氣相分離因子。依據拉烏爾定律,汽液相變的平衡常數可以由兩種不同氫同位素水在某一溫度下的飽和蒸汽壓值獲得。不同氫同位素水的飽和蒸汽壓有所不同,van Hook[9]基于凝聚狀態下蒸汽壓的同位素效應,采用計算方法獲得273~403 K下液態水的飽和蒸汽壓,并擬合成方程。根據飽和蒸汽壓比值獲得汽液相變的平衡常數。除了理論計算之外,加拿大和俄羅斯的研究人員同樣也對汽液相變過程的分離因子建立了經驗公式[3,7]。已知氣相反應和汽液相變兩步交換的平衡常數,最終也可以獲得水-氫交換總反應的分離因子,衡量不同同位素的分離效果。表1列出了HD/H2O、DT/D2O和HT/H2O三個體系常用的分離因子經驗公式。

為了對比經驗公式與理論計算值之間的差異,本工作對HD/H2O、DT/D2O和HT/H2O三體系的總分離因子αT進行比較,結果示于圖2。理論計算中反應的總分離因子是結合量化計算獲得的氣相反應分離因子和由飽和蒸汽壓比值獲得的汽液分離因子得到,從圖2可以看出:對于三種研究體系,兩經驗公式得到的較重氫同位素的總分離因子幾乎一致,并且經驗公式與理論計算得到的結果吻合較好,最大偏差僅為6.2%。因此,不管是理論計算還是經驗公式獲得的分離因子均能夠應用于模擬計算中。

2 結果與討論

2.1模型與方法驗證

基于平衡級模型,本工作首先依據模型建立計算方法,并驗證方法的可靠性,然后對HD/H2O、DT/D2O和HT/H2O體系氫同位素水-氫交換分離性能進行模擬研究。本工作以HD/H2O體系為例進行模型與方法的驗證。研究采用與文獻[1]中實驗一致的工藝條件為代表進行模擬驗證,工藝條件為:理論塔板數N=5;頂部進料H2O中D濃度為4×10-3(D原子百分數,下同);每一理論塔板壓降為133 Pa;氣液比G/L為1.0 ;交換柱底部進氣壓力為0.1 MPa和純H2(氣體中D濃度為0)。HD/H2O體系所采用的分離因子采用加拿大經驗公式(表1)。依據以上工藝條件進行模擬計算,獲得交換柱頂部氣相中D濃度隨溫度的變化,模擬結果示于圖3。由圖3可以看出,本工作建立模型獲得的模擬結果與文獻模擬結果一致,與實驗值偏差也較小。圖3中計算結果與實驗結果存在一定的差異。主要是因為計算值是在理想情況下獲得的,即氣液相完全接觸,并充分反應達到平衡。而實際條件下,尤其是在溫度比較高時,汽相增多會導致部分水汽在催化劑內部濃集,從而交換反應不能進行完全。除此之外,相比計算模擬的理想性,實際中反應的進行受到溫度梯度、質量傳遞速率等諸多不確定因素的影響。但是,所建立的模型與方法仍然是可靠的。

表1三種體系分離因子經驗公式

Table 1Empirical equations for the calculation of separation factors for three systems

體系國家分離因子HD/H2O加拿大[7]lnαg=-0.2735+449.2T+2380T2,lnαv=-0.0592-80.3T+25490T2,lnαT=-0.2143+368.8T+27870T2俄羅斯[3]lnαT=-0.1636+333.7T+33840T2日本[10]lnαg=-0.297+472.0THT/H2O加拿大[7]lgαg=0.292lgT+336.5T-1.055,lnαv=-0.00971-47.98T+23122T2,lnαT=-2.426+774T+0.292lnT俄羅斯[3]lnαT=-2.4264+718.2T+24589T2+0.292lnT日本[10]lnαg=-0.434+673.0TDT/D2O加拿大[7]lnαT=-0.1474+191.5T俄羅斯[3]lnαT=-0.1974+211.1T日本[10]lnαg=-0.132+170.0T

注:αg、αv和αT分別代表三體系中氣相反應、汽液相變和總反應中兩種氫同位素中較重氫同位素的分離因子

■——理論計算,□——加拿大經驗公式,△——俄羅斯經驗公式(a)——HD/H2O,(b)——DT/D2O,(c)——HT/H2O圖2 三種研究體系基于經驗公式與理論計算的分離因子對比圖Fig.2 Separation factors based on theoretical modelling and empirical equation for three systems

▲——本工作模擬值,●——文獻[1]模擬值,■——文獻[1]實驗值圖3 HD/H2O體系交換柱頂部氣相中D濃度隨溫度變化模擬值與實驗值對比Fig.3 Comparison of the simulated and experimental deuterium concentration at the column top for HD/H2O

2.2操作溫度的影響

本工作針對HD/H2O、DT/D2O和HT/H2O三個不同體系,研究了在293~353 K溫度范圍內,交換柱頂部氣相中較重氫同位素濃度隨溫度的變化。各體系在交換柱頂部氣相中較重氫同位素濃度分別代表:HD和H2混合氣中D濃度(HD/H2O);DT和D2混合氣中T濃度(DT/D2O);HT和H2混合氣中T濃度(HT/H2O)。在各研究體系中,交換柱底部進氣壓力均為0.1 MPa,進氣組成均為純H2(HD/H2O)、D2(DT/D2O)和H2(HT/H2O),每一理論塔板壓降為133 Pa;同樣頂部液體進料分別為HDO(D原子百分數為4×10-3)、DTO(T原子百分數為4×10-3)和HTO(T原子百分數為4×10-3)。HD/H2O 和HT/H2O體系所采用的分離因子來自于加拿大經驗公式;對于DT/D2O體系,沒有經驗公式獨立計算氣相反應和汽液相變的分離因子,因此DT/D2O體系分離因子來自于理論計算。模擬結果示于圖4。圖4結果表明:對于三種研究體系,交換柱頂部氣相中較重氫同位素濃度隨著溫度升高,均是先增加后降低的趨勢,存在最優的操作溫度343 K[11];在343 K時,HD中D(HD/H2O體系)和 DT中T(DT/D2O體系)的濃度高于HT中T(HT/H2O體系)的濃度,形成這種差別的主要原因是HT/H2O體系中氣相反應的分離因子比較大(4.63),高于HD/H2O體系(3.03)和DT/D2O體系(1.52),這取決于體系中分子熱力學性質的差別。進一步從統計熱力學角度分析,熱力學性質的差別取決于各分子的配分函數[8],分子配分函數之間的差別決定了熱力學性質的差別。

■——HD/H2O,●——DT/D2O,▲——HT/H2O圖4 三種研究體系交換柱頂部氣相中較重氫同位素濃度隨溫度變化Fig.4 Heavier hydrogen concentration at the column top for three systems as a function of temperature

2.3氣液比的影響

針對三種不同體系,接著研究了氣液比G/L從1.0增大到3.0的條件下,交換柱頂部氣相中較重氫同位素濃度在溫度與氣液比同時變化下的趨勢,結果示于圖5。

如圖5所示,在不同氣液比條件下,三體系中交換柱頂部氣相中較重氫同位素濃度均隨著溫度升高,呈現先增加后降低的趨勢,并且也存在一個最優的操作溫度。同時,最優操作溫度均隨著氣液比的增大呈現降低的趨勢,從343 K降低到323 K,在氣液比為1.5和2.0時,最優操作溫度均為333 K。

在氣液比從1.0增大到3.0的條件下,HD/H2O、DT/D2O中交換柱頂部氣相中較重氫同位素濃度分別下降45%和55%,而HT/H2O體系氣相中T的濃度降低36%,相比而言,降低程度較小。這說明在氣液比增大的過程中,HD/H2O、DT/D2O分離效果受溫度影響比較大,而HT/H2O體系的分離效果受溫度的影響較小一些。

G/L:■——1.0,●——1.5,▲——2.0,▼——3.0(a)——HD/H2O,(b)——DT/D2O,(c)——HT/H2O圖5 不同氣液比G/L條件下各研究體系交換柱頂部氣相中較重氫同位素濃度隨溫度變化Fig.5 Heavier hydrogen concentration at the column top for three systems as a function of temperature with various G/L ratios

2.4理論塔板數的影響

由以上分析可以得到對于三種體系,獲得最大分離效果的工藝條件為:操作溫度343 K和氣液比1.0。在此工藝條件下,進一步研究理論塔板數對分離效果的影響,交換柱頂部氣相中較重氫同位素濃度隨著理論塔板數N的變化示于圖6。從圖6可以看出,交換柱頂部氣相中較重氫同位素濃度均隨著理論塔板數增加先非線性升高,后逐漸趨于穩定。三種體系呈現的變化趨勢是一致的。HD/H2O和DT/D2O體系在理論塔板數為25的情況下趨于穩定而HT/H2O在理論塔板數為15的情況下就已經趨于穩定,這主要是因為HT/H2O體系的分離因子比較大,比HD/H2O和DT/D2O體系更快達到平衡。因此,在獲得最大的分離效果的要求下,HT/H2O體系需要的理論塔板數比HD/H2O和DT/D2O體系少。

■——HD/H2O,●——DT/D2O,▲——HT/H2O圖6 三種研究體系交換柱頂部氣相中較重氫同位素濃度隨理論塔板數的變化Fig.6 Heavier hydrogen concentration at the column top for three systems as a function of number of stages

2.5濃度分布

通過前述分析可以得到在獲得最大分離效果時,三種體系優化的工藝條件為:氣液比1.0,操作溫度343 K,理論塔板數25 (HD/H2O和DT/D2O)和15 (HT/H2O)。圖7給出了在優化的工藝條件下,各研究體系氣相中較重氫同位素濃度在交換柱內分布。從圖7可以看出,對于HD/H2O和HT/H2O體系,隨著理論塔板數的增加,即隨著交換柱高度的增加,各級理論塔板上的氣相中較重氫同位素濃度呈先非線性增加、后趨于穩定的變化規律。而對于DT/D2O體系,雖然也呈現非線性變化,但是濃度分布曲線的斜率隨著塔板數增加而增大。由前述方程(6)中可以得到,濃度分布曲線的斜率與系數AN成正比,三種體系的A值分別為:0.86(HD/H2O)、1.03(DT/D2O)和0.74(HT/H2O)。DT/D2O體系的A值大于1,因此對于DT/D2O體系,濃度分布曲線的斜率呈現增大的趨勢。系數A與體系的溫度、分離因子和液態水的飽和蒸汽壓有關。

■——HD/H2O,●——DT/D2O,▲——HT/H2O圖7 三種研究體系氣相中較重氫同位素濃度在交換柱內的分布Fig.7 Profiles of the heavier hydrogen concentration in the gas phase along the column

3 結 論

綜上所述,本工作采用數值模擬的方法,通過對HD/H2O、DT/D2O和HT/H2O三種氫同位素體系的分離性能對比研究,獲得結論如下:

(1) 交換柱頂部氣相中較重氫同位素濃度隨著溫度升高,均呈現先增加后降低的趨勢,存在一個最優的操作溫度343 K;

(2) 在氣液比從1.0增大到3.0的條件下,最優操作溫度均隨著氣液比的增大呈現降低的趨勢,并且HD/H2O、DT/D2O交換效果受溫度影響比較大,而HT/H2O體系的交換效果受溫度的影響相對較小;

(3) 在獲得最大分離效果的要求下,HT/H2O體系需要的理論塔板數比HD/H2O和DT/D2O體系少,同時,在優化工藝條件下,HD/H2O和HT/H2O體系氣相中較重氫同位素濃度隨著交換柱高度增大呈先非線性增加、后趨于穩定的相同變化規律,即濃度分布曲線的斜率呈現減小的趨勢,而DT/D2O體系,濃度分布曲線斜率呈現增大的趨勢。

[1]Ye Linsen, Luo Deli, Tang Tao, et al. Process simulation for hydrogen/deuterium exchange in a packed column[J]. International Journal of Hydrogen Energy, 2014, 39(12): 6604-6609.

[2]Kinoshita M. Simulation procedure for multistage water/hydrogen exchange column for heavy water enrichment accounting for catalytic and scrubbing efficiencies[J]. J Nucl Sci Technol, 1985, 22(5): 398-405.

[3]Fredorchenko O A, Alekseev I A, Trenin V D. Computer simulation of the water and hydrogen distillation and CECE process and its experimental verification[J]. Fusion Technol, 1995, 28 (2): 1485-1490.

[4]Kim K R, Ahn D H, Lee H S, et al. Rigorous analysis for design and performance of the water/hydrogen isotopic exchange column in a tritium removal facility[J]. J Ind Eng Chem, 1997, 3(4): 310-317.

[5]夏修龍.聯合電解催化交換系統HD/H2O和HT/H2O體系模擬[J].原子能科學技術,2006,40(Suppl):37-40.

[6]Isomura S, Kaetsu H, Nakane R. Deuterium separation by hydrogen-water exchange in multistage exchange column[J]. J Nucl Sci Technol, 1980, 17(4): 308-311.

[7]Rolston J H, den Hartog J, Bulter J P. The deuterium isotope separation factor between hydrogen and liquid water[J]. J Phys Chem, 1976, 80(10): 1064-1067.

[8]鐘正坤,張莉,孫穎,等.氫-水同位素交換分離因子理論計算[J].原子能科學與技術,2004,38(2):148-151.

[9]van Hook W A. Vapor pressures of the isotopic waters and ices[J]. J Phys Chem, 1968, 72(4): 1234-1244.

[10]Yamanishi T, Okuno K. A computer code simulation multistage chemical exchange column under wide range of operating conditions, JAERI-Data/Code 96-028[R]. Japan: Japan Atomic Energy Research Institute, 1996.

[11]Kim K R, Lee M S, Paek S, et al. Operational analysis of a liquid phase catalytic exchange column for a detritiation of heavy water[J]. Sep Purif Technol, 2007, 54(3): 410-414.

Simulation Study on H2O-H2Catalytic Exchange Capability for Hydrogen Isotopes Separation of HD/H2O, DT/D2O and HT/H2O

WU Dong, RUAN Hao, HU Shi-lin, YIN Yu-guo, ZHANG Li,HUANG Deng-gao, DOU Qin-cheng

China Institute of Atomic Energy, P. O. Box 275(83), Beijing 102413, China

Computational simulation, as a strong and valuable complement approach to experimental study, can play a significant role in hydrogen isotope separation nowadays. In this work, a systematic study was conducted on H2O-H2catalytic exchange capability for hydrogen isotopes separation of HD/H2O, DT/D2O and HT/H2O based on numerical modelling. The results indicate that the highest separation capability is obtained at 343 K, an optimum temperature for three systems. Additionally, the optimum temperatures decrease down to 323 K with increasing gas-liquid ratio up to 3.0, however, the HT/H2O system show less dependence on temperature considering the lower reduction of the separation ability. The number of theoretical stages for HT/H2O system is less than the other two ones as presenting best separation performance. The profile of hydrogen concentration along the column behaves differently for three systems under the technical conditions in this study.

H2O-H2exchange; hydrogen isotopes separation; numerical simulation; process conditions

2015-09-09;

2016-01-02

吳棟(1985—),男,山東聊城人,博士,從事氫同位素分離研究

O643.14

A

0253-9950(2016)04-0200-07

10.7538/hhx.2016.38.04.0200

猜你喜歡
理論體系研究
FMS與YBT相關性的實證研究
堅持理論創新
當代陜西(2022年5期)2022-04-19 12:10:18
神秘的混沌理論
遼代千人邑研究述論
理論創新 引領百年
構建體系,舉一反三
相關于撓理論的Baer模
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
“曲線運動”知識體系和方法指導
主站蜘蛛池模板: 亚洲中文字幕在线一区播放| 精品三级在线| 欧美日韩91| 久久精品国产一区二区小说| 国产一区在线观看无码| 久久人午夜亚洲精品无码区| 99久久婷婷国产综合精| 精品无码专区亚洲| 久久无码免费束人妻| 国产精品网拍在线| 一本一道波多野结衣av黑人在线| 国产99在线观看| 久久夜色精品国产嚕嚕亚洲av| 久久精品国产亚洲AV忘忧草18| 亚洲成人在线免费观看| 精品久久国产综合精麻豆| 亚洲国产欧美国产综合久久| 3p叠罗汉国产精品久久| 免费人成网站在线高清| 欧美日韩国产成人高清视频| 久久无码av一区二区三区| 国产迷奸在线看| 日韩欧美中文字幕在线精品| 乱系列中文字幕在线视频 | 久久久久国产精品熟女影院| 伊人久久综在合线亚洲91| 日韩欧美中文字幕在线韩免费 | 97青草最新免费精品视频| 日韩中文无码av超清| 国产精品刺激对白在线| 亚洲自拍另类| 国产精欧美一区二区三区| 91免费国产高清观看| 免费国产好深啊好涨好硬视频| 国产成人精品无码一区二| 亚洲一级毛片在线观播放| 欧美亚洲中文精品三区| 亚洲中字无码AV电影在线观看| 国产浮力第一页永久地址| 成人午夜网址| 久青草网站| 国产精品女人呻吟在线观看| 色哟哟色院91精品网站| 在线另类稀缺国产呦| a亚洲视频| 亚洲精品国产日韩无码AV永久免费网| 亚洲人成网址| 黄色成年视频| 91成人精品视频| 色网站免费在线观看| 国产成人亚洲精品色欲AV| 性欧美久久| 久久精品视频亚洲| 四虎永久免费在线| 2020极品精品国产| 亚洲性影院| 色屁屁一区二区三区视频国产| 熟妇丰满人妻av无码区| 国产在线一二三区| 日韩少妇激情一区二区| 日韩高清成人| 91亚洲国产视频| 久久久无码人妻精品无码| 香蕉eeww99国产在线观看| 国产视频资源在线观看| 美女高潮全身流白浆福利区| 亚洲美女久久| 国产地址二永久伊甸园| 色AV色 综合网站| 国产成人精品18| 国产主播一区二区三区| 97国产在线播放| 国产午夜不卡| 91娇喘视频| www.精品国产| 超薄丝袜足j国产在线视频| 精品黑人一区二区三区| 综合色婷婷| 中字无码精油按摩中出视频| 国产男女免费视频| 欧美日本中文| 亚洲中文字幕在线观看|