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

雙催化層固定床甲烷化反應器CFD模擬

2015-08-20 06:15:38趙靜張亞新冉文燊程源洪
化工學報 2015年9期
關鍵詞:催化劑模型

趙靜,張亞新,冉文燊,程源洪

(新疆大學化學化工學院,新疆 烏魯木齊 830046)

引 言

在煤制天然氣過程中,甲烷化技術是其關鍵技術之一。我國天然氣使用量占總能源的比例不及世界平均水平,隨著能源清潔要求越來越高,天然氣在能源結構中的地位日益突出。甲烷化反應作為放熱可逆的催化反應,在反應過程中溫度的分布對整個反應過程有極大的影響。目前,多數研究集中在催化劑的性能、催化劑載體、催化劑反應動力學性能、甲烷化工藝和反應器進出口參數之間的關系。

Wan 等[1]對Ru/MnO 催化劑在CO2/H2甲烷化反應中的催化特性做了全面研究。楊霞等[2-5]、張加贏等[6-8]對催化劑載體、助劑的特性進行了詳細分析。但小東等[9-12]使用ASPEN PLUS 軟件對甲烷化工藝進行了詳細分析,明確了出口參數與入口參數的關系。目前國內諸多學者對耐硫甲烷化反應催化劑進行了本征動力學分析[13-16],多選擇冪函數形式的經驗速率方程;對甲烷化反應器的數值模擬也局限于一維模型,多為平推流,不考慮徑向溫度變化[17-19]。實際上,甲烷化過程不僅與選擇的催化劑有關,還與許多因素相關。由于甲烷化反應是強放熱反應,對甲烷化反應器在反應狀態下內部流場的研究,特別是軸向、徑向溫度場特征的研究,對防止“飛溫”、保證甲烷產量有重要意義。

1 甲烷化反應CFX 模型的建立

1.1 甲烷化反應原理

一般情況甲烷化反應主要包括

反應過程中CO 對CO2的反應具有抑制作用[12],當CO 的量降到很小時,CO2反應才能有明顯的反應速率,在甲烷化反應體系中,從反應器入口到出口CO 的量逐漸減少,但始終存在,因此CO2的甲烷化反應量很少,本文只考慮式(1)。GCC 和MCR催化劑中反應速率方程式為

1.2 反應器模型及參數

圖1 為反應器內結構。反應器直徑為4900 mm,高度為9200 mm,采用外循環換熱,因此,反應器可視為絕熱。催化劑分兩層,GCC 催化劑層和MCR催化劑層,第1 層催化劑中反應較慢,為第2 層催化劑中大量反應做準備,氣體在GCC 催化劑中反應結束時溫度達到330℃,使反應氣溫度升高到MCR 催化劑的活性溫度范圍。兩種催化劑結合使用更有利于甲烷化反應。反應器下部剖面線部分是耐火材料做的襯里,熱導率小,有效避免了反應器溫度過高。本文采用三維擬均相多孔介質模型,對甲烷化反應器進行CFD 模擬。

圖1 絕熱甲烷化反應器結構Fig.1 Structure of adiabatic methanation reactor

表1 模型參數(從上到下)Table 1 Parameters of model (from top)

1.3 反應過程中的數學模型

1.3.1 多孔介質模型 床層可由多孔介質等效,由于催化劑顆粒并不均勻,造成孔隙率的不均勻,但孔隙率分布呈現一定的規律性。在壁面周圍孔隙率接近1,由壁面到中心逐漸減小到固定值。對于不同形狀的催化劑的孔隙率分布,許多學者提出了有關方程。Giese[20]研究得出孔隙率沿徑向的分布可以用以下函數來描述

孔隙率和分散系數都隨半徑不同而變化。熱量和質量平衡公式如下:

質量守恒

熱量平衡

1.3.2 k-ε模型 大多數CFD 研究選擇的湍流模型側重點各有不同,本文采用目前應用最為廣泛的標準k-ε速度與長度兩方程模型,其中k 為湍動能,表示速度波動的變化量;ε是湍動能耗散,表示速度波動耗散的速率。ε方程為

對于不可壓縮流體,不計重力影響,其中k 和ε方程分別定義為

1.3.3 熱能模型 忽略由于黏度引起的內部熱量,通過對動量方程的積分得到熱能方程

1.4 網格劃分及邊界條件

網格劃分工具采用ICEM CFD專用的流體網格劃分軟件,流體區域采用四面體網格。由于四面體網格不能很好地將圓弧特征描述出來,在邊緣部分使用三棱柱網格,既能很好表征圓弧特征,又能減少網格數量。設定入口質量流量、入口氣體組成、入口氣體溫度、出口壓力,壁面無滑移,孔隙率隨半徑變化,設置化學反應,對模型進行穩態模擬。另外,催化劑失活溫度是700℃,本文中假設溫度達到700℃或700℃以上催化層的反應速率為0。

2 計算結果及分析

2.1 模擬結果特征參數對比

對模擬結果與設計數據進行對比分析,結果見表2。

表2 模擬結果與設計數據對比Table 2 Simulation result compared with design data

通過以上對比,各主要參數模擬值和設計值的誤差均≤5%,可以認為模擬方法的正確性和結果的可信性。

2.2 模擬結果及分析

反應器入口氣體組分是CO 28%、H26.41%、CH434.5%、H2O 25.77%、N24.32%。圖2 為沿中心縱截面的溫度分布模擬結果。由圖可明顯看出,在GCC 催化層中沿徑向溫度分布較均勻,中心與近壁面溫差較小,隨著反應加劇,在MCR 催化層中沿徑向溫度分布不均勻性更明顯,催化劑出口處由于形狀突變和下游流體域不一致引起的溫度不均勻性尤其突出。

圖3 為不同橫截面上溫度云圖,圖4 為圖3 對應截面上溫度沿徑向的分布曲線,截面取離原點分別為1、1.9、2.9、3.7、4.5 m 依次沿床層由上向下。由圖4 可以看出,距原點越遠,對應截面上徑向溫度分布越不均勻。1.9 m 處截面處于GCC 催化劑出口,在此截面上溫度沿徑向波動很小,幾乎均勻分布。只有當截面接近MCR 催化劑出口時,截面徑向溫度才有所差別,特征為溫度呈軸對稱分布,由中心向壁面逐漸升高。

圖2 沿中心軸截面溫度分布Fig.2 Temperature distribution along central axis

圖3 不同橫截面溫度云圖Fig.3 Temperature nephogram on different cross section

圖4 溫度沿徑向的分布Fig.4 Temperature distribution along radial curve

圖5 流速沿徑向的分布Fig.5 Velocity distribution along radial curve

圖5 為氣體由原點向下1、1.9、2.7、3.7 m 截面上流速模擬分析曲線。由圖5 看出,在催化劑床層出口附近,流速呈中心大、兩邊小趨勢,這是造成催化劑床層出口截面溫度徑向分布呈中心低、兩邊高現象的主要原因。在催化劑中心處反應集中,產熱多,同時中心氣體流速大,熱量交換快,帶走的熱量多,因此在催化劑床層出口處存在中心溫度低,近壁面溫度高的現象。

離原點1.9 m 處為GCC 催化劑層出口處,此處下游反應器橫截面縮小。此截面流速分布特點是:半徑約2 m 區域內流速基本均勻,半徑從2 m 到壁面間流速逐漸變小,此截面上的溫度分布也存在相應的不均勻現象。

3 影響因素分析

3.1 床層結構影響

通過對模擬結果的分析可知,床層直徑發生變化處對鄰近區域,尤其是上游區域影響較大。床層溫度越均勻,催化劑利用率越高,生成甲烷量越多。希望改變床層結構提高床層溫度均勻性,避免局部過熱出現“飛溫”現象。

3.1.1 床層等徑處理 圖6 為等徑模型,通過改變反應器直徑,使反應器變徑部位減少。等徑模型下反應器內溫度分布如圖7 所示,可以看出,等徑結構下溫度分布更均勻,更接近帶狀分布,只有在催化劑床層出口附近出現明顯的近壁面溫度升高現象。但是,等徑結構由于減少了變徑處氣體沿變徑邊界層的渦流擾動,造成出口甲烷含量(43.85%)低于原模型。

圖6 等徑模型Fig.6 Equal radius model

圖7 溫度分布Fig.7 Temperature profile with equal radius model

3.1.2 拆分催化層 在MCR 催化劑床層中間加支撐層,將MCR 催化層分成0.7 和1 m 兩段,加入的支撐層厚100 mm,拆分后的模型見圖8。對圖8中模型進行模擬,溫度場分布如圖9 所示。模擬結果顯示,溫度分布呈中心低、近壁面高的現象在第2 段催化劑中更加明顯。

圖8 MCR 段床層拆分處理Fig.8 MCR section layer separated model

圖9 拆分床層下溫度分布Fig.9 Temperature profile of MCR section layer separated model

起初希望在出現溫度分布不均勻的位置增加支撐分隔以使氣體得以緩沖、混合更均勻,有利于緩解催化劑出口處溫度分布不均現象。模擬結果顯示,這種拆分催化床層100 mm 間隔的方法對溫度分布均勻幾乎沒有效果,這可能是由于分段的支撐層位置不在溫度分布最不均勻位置或者拆分段間距太小,有必要對拆分間隔與溫度均勻化問題進一步分析。

3.1.3 延長催化層出口段支撐 鑒于床層底部催化劑出口處圓弧處理后對整個床層內各場分布的有利影響,對催化劑出口段支撐等徑部分進行延長(至500 mm),建立的模型如圖10 所示。對圖10 模型進行溫度場分析,得到溫度分布云圖如圖11 所示。可以看出,最高溫度較原模型的973℃下降到955℃,且溫度分布更加均勻,出口甲烷的平均質量分數由45.02%升高到45.21%,說明本結構有利于溫 度均勻分布和甲烷質量分數的提高。

圖10 催化劑出口段支撐延長處理Fig.10 End support extend model

圖11 支撐延長處理后溫度分布Fig.11 Temperature profile with end support extend model

3.1.4 改變床層結構的對比評價 本文共進行了MCR 催化劑等直徑處理、定催化劑體積下直徑為原徑1.2 倍和0.8 倍、MCR 催化層拆分分隔100 mm、催化層出口段支撐延長等反應器和床層結構改進,其中各結構對應的出口甲烷質量分數見表3。

各種模型下MCR 床層出口截面溫度徑向分布如圖12 所示。從圖可以看出,MCR 催化層的分層處理模型與原模型的催化層的出口處溫度分布幾乎沒有差異。直徑變為原來的1.2 倍時,溫度分布最不均勻,直徑為原0.8 倍時,模型中能量損耗大,一段溫升大,不利于設備安全。支撐段延長并做圓弧處理使催化劑層溫度分布最均勻,橫截面溫差最小,且出口甲烷質量分數有所提高,整個溫度場幾乎不存在超溫現象。等徑處理雖然降低了溫度,但是反應氣體的停留時間也變短,氣體反應時間變短,造成甲烷的產量減低。

表3 各結構對應的甲烷質量分數Table 3 Mass fraction of CH4 with different structure

圖12 結構改變對應的催化層出口溫度分布Fig.12 Outlet temperature profile corresponding to structure change

綜合考慮甲烷的產量、能耗、設備安全性,本文認為最優模型應為出口段延長并做圓弧處理的模型。

3.2 入口參數影響

影響甲烷化反應產物轉化率的因素眾多,從甲烷生成的速率方程看,影響甲烷轉化率的因素有溫度和反應物的分壓,溫度與反應速率呈指數關系,對反應速率的影響更大。此外,入口氣體的流量、入口壓力都是操作控制參數。因此,本文對入口溫度、入口流量和入口壓力進行研究,提出其相應的允許波動范圍。波動約束條件:反應器的工作溫度不能超過設計溫度400℃,床層內溫度不能超過催化劑的活性溫度范圍,催化劑失活率不超過5%,甲烷的產率不低于45%。

3.2.1 入口溫度影響 入口溫度在設計溫度左右取6 組數據,得到相應的床層最高溫度、催化劑失活率和出口甲烷的質量分數如表3 所示。

表中床層溫度隨入口溫度的升高而升高,甲烷的質量分數也略有升高,但是反應器上部壁面溫度可能超過設備設計溫度。

表4 各參數隨入口溫度的變化Table 4 Parameters changing with inlet temperature

入口溫度的上升對甲烷產率的影響是“雙刃劍”。一方面,入口溫度的升高必然使氣體進入催化劑層的溫度升高,由反應動力學可知,反應速率也隨著加大,生成的甲烷量也隨著增加;另一方面,隨著入口溫度不斷升高,催化劑失活率也開始提高,失活導致甲烷產率下降。從結果上看,甲烷產量增加,說明入口溫度帶來的反應加速大于催化劑失活造成的產率下降,提高進口溫度對提升產率起著主導作用。

入口溫度的升高,不僅使催化劑失活率升高,而且一段催化劑壁面溫度也升高,不利于設備安全。權衡甲烷的產率與催化劑失活率,甲烷質量分數不低于45%,催化劑失活率不高于5%,當入口溫度為262℃時,催化劑失活率超過5%,本文建議入口溫度應嚴格控制在253~262℃。

3.2.2 空速影響 設計參數下反應器的入口空速為9424.24 kg·m-3·h-1。本文取空速7500、8500、9000、10500 和11500 kg·m-3·h-15 組數據進行空速影響研究。

圖13 是在不同空速下,床層內出現最高溫度截面上提取的沿徑向的溫度分布曲線,由于床層在不同空速下最高溫度出現的截面不同,所以,提取的徑向溫度分布并不在床層同一截面上,與空速對應的不同截面截取位置依次是3.2、3.4、3.6、3.75、3.85、3.85 m,此處距離均為截面到床層頂部原點的距離。在所取截面上,溫度分布最不均勻,而且同一半徑上溫度也最高。空速越小溫度分布越不均勻,出口溫度越高,相應的甲烷的質量分數越高,且空速大時壓降增大,能耗增加。這種現象的理論解釋為:空速增大即入口質量流量增大,通過催化劑的氣體速度增大,原料氣在催化劑中的停留時間縮短,使甲烷的質量分數降低。流速增大氣體交換加快,反應生成的熱被氣體帶走,整個床層的溫度降低。

空速應穩定在一定的范圍內。當空速小于8500 kg·m-3·h-1時,催化劑會存在大范圍(>13.3%)失活,且第1 段催化劑溫度過高,引起設備安全問題;空速過大引起甲烷的質量分數過低,無法達到生產要求。所以建議空速應在9529~9345 kg·m-3·h-1之間,以保證催化劑活性和甲烷的產量。

圖13 反應器最高溫度截面處溫度隨空速的變化Fig.13 Temperature profile with changing space velocity in Tmax section

3.2.3 入口壓力影響 在實際操作中通過調節入口壓力來調節影響反應器的工藝參數,入口壓力和反應器內壓力損失決定反應器出口壓力,在CFX 設置中通過改變出口壓力反求入口壓力,保持其他條件不變,研究反應器出口參數等隨反應器入口壓力變化的規律。

催化劑床層出口溫度隨入口壓力變化曲線如圖14 所示。圖中看出,催化劑出口溫度隨著入口壓力的增大而升高,壓力越小溫度分布越均勻,最高溫度也越低,甲烷的出口質量分數越低;入口壓力越高,溫度分布越不均勻,出口處甲烷的平均質量分數越大。壓力升高單位體積氣體量增加,氣體的濃度增加,氣體的反應速率增大,床層溫度升高。

圖14 催化劑出口溫度隨入口壓力變化曲線Fig.14 Outlet temperature of catalyst changing with inlet pressure curve

入口壓力升高到3.273 MPa 時,催化劑失效體 積分數達6.08%,而且壓力升高對設備的要求升高,不利于設備安全和生產穩定。綜上所述,在其他條件保持設計值以內情況下,入口壓力應維持在3.120~3.251 MPa 范圍內。

4 結 論

本文通過對某工藝中雙段床層甲烷化反應器的CFD 模擬計算與分析,得到以下結論。

(1)建立了一種以數值模擬手段對反應狀態下多催化床層進行特征場研究的方法。

(2)通過ANSYS-CFX 軟件對甲烷化反應器進行傳質、傳熱和化學反應的模擬計算,得到甲烷化反應器內在反應狀態下溫度場、壓力場、速度場分布規律,為后續對甲烷化反應器催化床層的結構優化提供依據。

(3)對幾種不同床層結構進行分析,以結構安全為前提,得到了床層結構與溫度分布、出口甲烷產量的關系;確定了既有利于設備安全,也有利于甲烷產量的提高的相對優化的床層結構模型。

(4)通過對影響甲烷化反應速率的因素研究,明確入口參數與溫度分布以及甲烷產量之間的關系;提出了針對本工藝的反應器入口參數的允許波動范圍。

符 號 說 明

C1ε, C2ε——線性和平方阻力系數

dp——小球直徑,mm

k ——氣體湍動能,m2·s-2

pCO, pH2——分別為CO、H2的氣體分壓,Pa

R ——氣體常數,J·mol-1·K-1

r ——催化劑半徑,mm

rCH4——MCR 催化劑的反應速率,kg·s-1·kg-1

r′CH4——GCC 催化劑的反應速率,kg·s-1·kg-1

rtube——反應器半徑,mm

T ——氣體熱力學溫度,K

ε ——氣體湍動能耗散,m2·s-3

μ ——氣體動力黏度,Pa·s

ρ ——氣體密度,kg·m-3

ρbed——床層密度,kg·m-3

[1]Abu Bakar W A W, Toemen S, Ali R, Abd Rahim H F.Elucidation of active species over Ru/MnO catalyst on CO2/H2methanation reaction [J].Advances in Materials Physics and Chemistry, 2013, 3(2): 161-167.

[2]Yang Xia(楊霞), Tian Dayong(田大勇), Sun Shouli(孫守理), Sun Qi(孫琦).Effect of CeO2on the performance of nickel-based catalysts for methanation [J].Industry Catalysis(工業催化), 2014, 22(2): 137-143.

[3]Yang Xia(楊霞), Zheng Wentao(鄭文濤), Wang Guogao(汪國高), Sun Shouli(孫守理), Sun Qi(孫琦).Effect of MgO on catalytic performance of Ni/Al2O3catalyst for CO methanation [J].Modern Chemical Industry (現代化工), 2014, 34(1): 90-94.

[4]Yang Xia(楊霞), Tian Dayong(田大勇), Sun Shouli(孫守理), Sun Qi(孫琦).Influence of zirconia-alumina composite on catalytic performance of nickel-based catalysts for methanation [J].Chemical Industry and Engineering Progress(化工進展), 2014, 33(3): 673-678.

[5]Zhang Guoquan, Peng Jiaxi, Sun Tianjun, Wang Shudong.Effects of the oxidation extent of the SiC surface on the performance of Ni/SiC methanation catalysts [J].Chinese Journal of Catalysis, 2013, 34(3): 1745-1755.

[6]Zhang Jiaying(張加贏), Xin Zhong(辛忠), Meng Xin(孟鑫), Tao Miao(陶淼).Activity and stability of nickel based MCM-41 methanation catalysts for production of synthetic natural gas [J].CIESC Journal(化工學報), 2014, 65(1): 160-168.

[7]Zhang Xintang(張新堂), Yang Ying(楊影), Li Zhenxing(李振興), Sun Miaoyuan(孫淼元).The influence of composite support on the methanation activity of nickel-based catalysts [J].Journal of Shandong University of Science and Technology(山東科技大學學報), 2014, 33(1): 53-58.

[8]Li Xia(李霞), Yang Xiazhen(楊霞珍), Tang Haodong(唐浩東), Liu Huazhang(劉化章).Effect of supports on catalytic performance of nickel-based catalyst for methanation [J].Chinese Journal of Catalysis(催化學報), 2011, 32(8): 1400-1404.

[9]Dan Xiaodong(但小東), Yan Li(顏麗), Yang Yang(楊洋).Simulation and analysis of coke oven gas methanation process based on ASPEN PLUS [J].Technology and Development of Chemical Industry(化工技術與開發), 2014, 43(5): 52-57.

[10]He Yifu (何一夫).Simulation of methanation process based on ASPEN PLUS [J].Modern Chemical Industry(現代化工), 2012, 32(4): 107-109.

[11]Tian Liang(田亮), Jiang Da(蔣達), Qian Feng(錢鋒).Simulation and optimization of acetylene converter with decreasing catalyst activity [J].CIESC Journal (化工學報), 2012, 63(1): 185-192 .

[12]Bukur D B, Pan Zhendong, Ma Wenping, Jacobs G, Davis B H.Effect of CO conversion on the product distribution of a Co/Al2O3Fischer-Tropsch synthesis catalyst using a fixed bed reactor [J].Catal.Lett., 2012, 142(11): 1382-1387.

[13]Yu Jianguo(于建國), Yu Zunhong(于遵宏), Sun Xingyuan(孫杏元), Pan Huiqin(潘慧琴), Yu Guangsuo(于廣鎖).Macrokinetics of sulfur-tolerant methanation SDM-1 catalyst [J].Journal of Chemical Industry and Engineering(China) (化工學報), 1994, 45(1): 120-124.

[14]Yu Guangsuo(于廣鎖), Yu Jianguo(于建國), Sun Xingyuan(孫杏元), Pan Huiqin(潘慧琴), Yu Zunhong(于遵宏).Intrinsic kinetic study of KD306-type sulfur tolerant methanation catalyst [J].ECUST Journal (華東理工大學學報), 1999, 25(2): 15-17.

[15]Yu Jianguo(于建國), Lei Hao(雷浩), Wang Fuchen (王輔臣), Pan Huiqin(潘慧琴), Sun Xingyuan(孫杏元).Kinetic studies of methanation on SDM-1 sulfur-tolerant catalyst model for intrinsic kinetics [J].Journal of Fuel Chemistry and Technology (燃料化學學報), 1994, 22(3): 225-231.

[16]Yu Jianguo(于建國), Shen Caida(沈才大), Yu Guangsuo(于廣鎖), Sun Xingyuan(孫杏元), Pan Huiqin(潘慧琴).Kinetic studies of methanation on SDM-1 sulfur-tolerant catalyst model of global kinetics [J].Journal of Fuel Chemistry and Technology (燃料化學學報), 1994, 22(3): 232-238.

[17]Er-rbib H, Bouallou C.Methanation catalytic reactor [J].Comptes Rendus Chimie, 2014, 17(4): 701-706.

[18]Nijemeisland M, Dixon A G.Comparison of CFD simulations to experiment for convective heat transfer in a gas-solid fixed bed [J].Chemical Engineering Journal, 2001, 82(1): 231-246.

[19]Xu Jian, Wei Weisheng, Tian Aizhen, Fan Yu, Bao Xiaojun, Yu Changchun.Temperature profile in a two-stage fixed bed reactor for catalytic partial oxidation of methane to syngas [J].Catalysis Today, 2010, 149(4): 191-195.

[20]Giese M.Str?mung in por?sen Medien unter Berücksichtigung effektiver Viskosit?ten [D].Bavaria: TU München, 1997.

猜你喜歡
催化劑模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
直接轉化CO2和H2為甲醇的新催化劑
鋁鎳加氫催化劑在BDO裝置運行周期的探討
3D打印中的模型分割與打包
新型釩基催化劑催化降解氣相二噁英
掌握情欲催化劑
Coco薇(2016年2期)2016-03-22 02:45:06
V2O5-WO3/TiO2脫硝催化劑回收研究進展
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲性影院| 色婷婷狠狠干| 成人在线不卡| 国产亚洲精品自在久久不卡| 欧洲日本亚洲中文字幕| 伊人大杳蕉中文无码| 国产拍揄自揄精品视频网站| 欧美一级在线播放| 亚洲天堂视频在线观看免费| 国产9191精品免费观看| 精品少妇人妻无码久久| 青青草原国产精品啪啪视频| 国产99视频在线| 99re这里只有国产中文精品国产精品 | 国产欧美精品一区二区| 国产综合无码一区二区色蜜蜜| 波多野结衣一区二区三区AV| 亚洲第七页| 婷婷六月综合网| 九九久久精品免费观看| 97影院午夜在线观看视频| 欧美激情第一区| 国产91无码福利在线 | 国内老司机精品视频在线播出| 婷婷亚洲综合五月天在线| 久久久精品无码一二三区| 热re99久久精品国99热| 午夜激情福利视频| 国产一区三区二区中文在线| 亚洲国产成人麻豆精品| 午夜欧美理论2019理论| 在线欧美日韩国产| 亚洲免费播放| 欧美福利在线观看| 亚洲制服丝袜第一页| 中日无码在线观看| 国产全黄a一级毛片| 欧美成人亚洲综合精品欧美激情| 97在线视频免费观看| 激情国产精品一区| 久草视频中文| 97色婷婷成人综合在线观看| 欧美亚洲第一页| 久久96热在精品国产高清| 日韩AV无码一区| 亚洲福利一区二区三区| 污网站在线观看视频| 欧美色99| 在线欧美一区| 久久一级电影| 亚洲欧洲日本在线| 无码中文字幕精品推荐| 成年A级毛片| 亚洲精品国产精品乱码不卞 | 色综合婷婷| 国产激情无码一区二区三区免费| 国产在线一区二区视频| 777午夜精品电影免费看| 无码高潮喷水专区久久| 99视频在线免费观看| 亚洲午夜国产片在线观看| 国产精品美人久久久久久AV| 国产爽歪歪免费视频在线观看 | 亚洲成人黄色网址| 国产欧美精品一区aⅴ影院| 国内黄色精品| 二级特黄绝大片免费视频大片| 国产午夜不卡| 狠狠v日韩v欧美v| 国产精品刺激对白在线| 亚洲中文字幕无码mv| 亚洲AⅤ波多系列中文字幕| 美女一区二区在线观看| 伊伊人成亚洲综合人网7777| 亚洲欧美日韩久久精品| 欧美成人A视频| 国产日韩丝袜一二三区| 在线亚洲精品自拍| 人妻丰满熟妇αv无码| 高清国产va日韩亚洲免费午夜电影| 国产成人综合日韩精品无码首页| 中文字幕欧美成人免费|