亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        熱等離子體裂解丙烷制乙炔數(shù)學模擬

        2016-09-26 08:09:13蘇寶根吳劍驊聞光東房建威邢華斌任其龍
        化學反應工程與工藝 2016年4期
        關鍵詞:丙烷等離子體反應器

        蘇寶根,吳劍驊,聞光東,馬 杰,房建威,邢華斌,任其龍

        浙江大學生物質(zhì)化工教育部重點實驗室 化學工程與生物工程學院,浙江 杭州 310027

        熱等離子體裂解丙烷制乙炔數(shù)學模擬

        蘇寶根,吳劍驊,聞光東,馬 杰,房建威,邢華斌,任其龍

        浙江大學生物質(zhì)化工教育部重點實驗室 化學工程與生物工程學院,浙江 杭州 310027

        為了對反應器的停留時間進行合理優(yōu)化,將旋轉(zhuǎn)弧等離子體反應器視為一維平推流反應器網(wǎng)絡模型,結(jié)合裂解反應動力學模型與反應器流動模型,采用 CHEMKIN-PRO對丙烷的裂解過程進行數(shù)值模擬,用于分析熱等離子體反應器內(nèi)丙烷的裂解過程中產(chǎn)物的濃度分布及溫度分布情況。反應動力學模型分別采用均相反應動力學模型和非均相反應動力學模型。模擬結(jié)果表明,包含結(jié)焦模型的非均相反應動力學模型與實驗結(jié)果表現(xiàn)出更好的一致性,隨著反應器長度的增加,乙炔濃度存在最佳點。通過降低反應器的停留時間至1.0 ms以下,能有效提升C2H2收率。

        數(shù)學模擬 熱等離子體 丙烷 裂解 乙炔

        熱等離子體具有高溫、高焓以及富含活性粒子等特性,能夠在毫秒級的時間內(nèi)完成反應,實現(xiàn)物質(zhì)的高效轉(zhuǎn)化。熱等離子體技術已被廣泛運用于煤裂解制乙炔[1]、煤氣化[2]、碳烷烴轉(zhuǎn)化[3]、二氧化碳重整[4]以及生物質(zhì)氣化[5]等研究中。因不產(chǎn)生CO2及耗水量少,熱等離子體裂解煤及低碳烷烴制乙炔工藝被視為能代替?zhèn)鹘y(tǒng)電石法的既環(huán)保又高效的新工藝。

        由于熱等離子體裂解制乙炔的反應過程在毫秒級時間內(nèi)完成,產(chǎn)物需要快速淬冷以防止乙炔的進一步轉(zhuǎn)化,因此,對該反應的動力學研究對于熱等離子體反應器的設計及優(yōu)化顯得尤為重要。然而由于等離子體裂解在瞬時完成反應,對其動力學的實驗研究較難實現(xiàn)?,F(xiàn)階段相關的熱等離子體裂解動力學研究主要集中于甲烷裂解制乙炔[6,7]及煤的裂解過程[8]。熱等離子體裂解體的熱力學平衡研究表明,不包含固態(tài)碳的均相熱力學平衡計算結(jié)果與裂解的實驗結(jié)果更為相符[9]。研究者們普遍認為在熱等離子體反應器內(nèi)的停留時間為毫秒級,這使得反應器內(nèi)氣固相平衡難以達到,整個反應體系可認為是單一氣相平衡狀態(tài)。故有研究者在動力學研究中采用均相反應動力學模型進行模擬[6,7]。Slovetskii等[6]選擇了29種物質(zhì)及74個氣相反應用于描述裂解過程中甲烷的轉(zhuǎn)化及主要產(chǎn)物(C2H2及H2)的形成,并結(jié)合反應器內(nèi)的傳質(zhì)及傳熱過程建立了一個甲烷裂解動力學的數(shù)學模型。模擬結(jié)果表明,產(chǎn)物中甲烷濃度的計算值與實驗值存在較大偏差,認為主要誤差來源于等離子體炬用于加熱氫氣的實際功率、氫氣與甲烷的氣體流量、反應時間及最初混合物中的氫碳比等。由于在實際裂解反應中,焦的生成無法避免,因此有研究者采用非均相反應動力學模型進行模擬。Holmen等[10]提出了一個簡化的甲烷等離子體裂解的表觀動力學模型,并利用該簡化模型對甲烷的裂解過程進行計算,得到了反應體系中各組分隨反應時間的分布情況。Fincke等[11]建立了包含29種物質(zhì),87個氣相反應及2個結(jié)焦反應的動力學模型,用于研究甲烷的裂解反應動力學。在此基礎上,利用CHEMKIN 3.6中一維平推流反應器模擬了不同條件下甲烷裂解結(jié)果。利用該方法計算的模擬結(jié)果與實驗值的偏差較小。

        由于在等離子體環(huán)境中丙烷在熱力學上并不能穩(wěn)定存在,因此早期的動力學研究模型中并不包含丙烷,無法直接用于丙烷裂解過程的動力學模擬研究。為此,本工作用一維平推流反應器網(wǎng)絡模型對MW級旋轉(zhuǎn)弧等離子體反應器進行簡化處理,分別結(jié)合均相反應動力學模型及非均相反應動力學模型(包含結(jié)焦模型)進行計算模擬。均相反應數(shù)據(jù)庫包含198種物質(zhì)與762個氣相基元反應,能適用于多種化合物的裂解過程研究。結(jié)焦模型則采用了燃燒過程中較為公認的多環(huán)芳烴成核,氫消除乙炔加成(HACA)[12,13]及多環(huán)芳烴縮聚[14]的固體焦表面生長機理。通過對比丙烷裂解的模擬結(jié)果與實驗值,篩選出更加合理的動力學模型,用于指導工藝條件和反應器結(jié)構(gòu)優(yōu)化。

        1 丙烷等離子體裂解過程模型的建立及描述

        1.1反應器模型

        旋轉(zhuǎn)弧等離子體反應器的主要結(jié)構(gòu)如圖1所示,它由棒狀陰極、套筒狀陽極、淬冷段及勵磁線圈組成。反應器內(nèi)電弧在外加磁場的作用下高速旋轉(zhuǎn),形成均勻的等離子體。丙烷與氫氣從旋轉(zhuǎn)電弧上方通過,直接與等離子體充分接觸,快速完成裂解反應,并在淬冷器的作用下迅速冷卻至100℃以下,得到穩(wěn)定的裂解氣。

        圖1 等離子體炬反應器Fig.1 Schematic of rotating arc plasma reactor

        整個反應器主要為管式結(jié)構(gòu),采用CHEMKIN-PRO中的平推流反應器模塊,建立簡化的反應器網(wǎng)絡,用于模擬丙烷的等離子體裂解過程,簡化的反應器網(wǎng)絡模型如圖2所示。

        圖2 CHEMKIN-PRO的反應器網(wǎng)絡模型Fig.2 Reactor network model in CHEMKIN-PRO

        模型建立過程中的基本假設如下:1)反應器內(nèi)流體的流動過程被視為一維穩(wěn)態(tài)平推流,所有的物質(zhì)參數(shù)在徑向方向上均勻分布;2)等離子體的初始溫度通過氣體的比焓(輸入功率/氣體流量)換算得到;3)實際淬冷過程中淬冷水與產(chǎn)品氣之間的反應過程被忽略,淬冷過程理想化。利用該模型研究不同功率條件下丙烷的等離子體裂解過程,模擬所用的實驗條件如表 1所示。表中,H2,C3H8和Ar的體積百分數(shù)分別為99.50%,94.00%和99.99%,H2,C3H8和Ar的流量均為標準狀態(tài)下的流量。

        表1 丙烷裂解條件Table 1 Conditions of propane pyrolysis

        1.2動力學模型

        均相反應動力學模型中包含有198種物質(zhì)與762個氣相基元反應,用于描述包括丙烷裂解、乙炔生成以及作為結(jié)焦前體多環(huán)芳烴的生成過程。其中C1~C6的基元反應主要來源于JetSurF V2.0模型[15]與USC Mech V2.0模型[16]。USC Mech V2.0模型描述了C1~C4烴類的高溫裂解與氧化過程,JetSurF V2.0模型進一步完善了包括 C5~C12中直連烷烴、環(huán)烷烴以及單烷基取代的環(huán)烷烴在高溫下的裂解與氧化過程。多環(huán)芳烴(PAHs)的生長過程與聚炔烴(C2nH2)的生長過程則參照氫消除乙炔加成(HACA)中所涉及的基元反應進行設置。CHEMKIN-PRO中反應速率常數(shù)的計算基于修正的阿累尼烏斯方程,如下式所示:

        式中:A為指前因子,cm3/(mol.s);T為反應溫度,K;n為溫度的指數(shù)修正項;Ea為活化能,kJ/mol;R為摩爾氣體常量,8.314 J/(mol.K)。

        非均相反應動力學模型在均相反應動力學模型的基礎上,引入了固體焦的形成模型。焦的形成模型通常分為兩個過程,包括焦的成核過程與表面生長過程。其中,多環(huán)芳烴[17,18]與聚炔烴[19]在燃燒及裂解的研究中被廣泛視為結(jié)焦前體。本工作選擇芘(A4)[20,21]作為結(jié)焦前體,焦的成核速率可視為芘的分解速率。焦的表面生長過程包含了氫消除乙炔加成以及多環(huán)芳烴縮聚反應。固體焦形成反應如表 2所示。在固體焦表面的活性位點(α)對焦的生長速率有較大的影響,本工作選擇Wen等[22]報道的α值(為1.66×10-7mol/cm2)用于固體焦生長速率的計算。

        表2 固體焦的形成反應Table 2 Reaction of soot formation

        2 結(jié)果與討論

        2.1均相反應動力學模型模擬結(jié)果分析

        采用均相反應動力學模型對丙烷的等離子體裂解過程進行模擬,實驗條件如表1所示。不同實驗條件下的模擬計算結(jié)果與實驗值列于表3,表中模擬I為采用均相模型的模擬結(jié)果。反應器內(nèi)的溫度分布及主要產(chǎn)物CH4,C2H4和C2H2的濃度分布如圖3所示。

        圖3 均相反應動力學模型模擬結(jié)果中CH4,C2H4和C2H2的濃度分布與溫度分布Fig.3 Distributions of CH4, C2H4and C2H2concentrations and temperature from the simulation of homogeneous reaction kinetic model

        由圖可以看出,反應器內(nèi)的溫度由于裂解反應而迅速下降,并在冷卻水的影響下進一步降低,不同實驗條件下所得的淬冷前溫度分別為1 695,1 795及2 371 K。在反應段前端約為0.02~0.10mm處,CH4與C2H4的濃度達到最大值,隨后兩者的濃度急劇下降至平衡值。C2H2濃度在反應段前端急劇上升,在10~20mm處,C2H2濃度趨于穩(wěn)定,而隨著反應段長度的增加,C2H2濃度的變化幅度不明顯。對比不同功率條件下的模擬結(jié)果,隨著輸入功率的增加,C2H2濃度在反應段內(nèi)的分布趨勢發(fā)生改變。當輸入功率為614.9和692.8 kW時,C2H2濃度在反應段內(nèi)達到平衡后,隨著反應器長度的進一步增加,C2H2濃度有所下降;當輸入功率增加至794.2 kW 時,C2H2濃度在反應段內(nèi)達到平衡后,隨著反應器長度增加C2H2濃度一直上升。在淬冷段中,C2H2的濃度變化不大。由表3可以看出,反應器出口的最終結(jié)果顯示,C2H2產(chǎn)率隨輸入功率的增加而增加,然而實驗結(jié)果顯示C2H2產(chǎn)率存在最佳值。均相反應動力學模型模擬結(jié)果所得的C2H2濃度及產(chǎn)率均高于實際值??梢娪镁喾磻獎恿W模型進行模擬,與實際的反應過程存在較大的偏差。

        2.2非均相反應動力學模型模擬結(jié)果分析

        采用非均相反應動力學模型對丙烷的等離子體裂解過程進行模擬,實驗條件如表 1所示。不同實驗條件下的模擬計算結(jié)果與實驗值列于表3,表中模擬II為采用非均相模型的模擬結(jié)果。反應器內(nèi)的溫度分布及主要產(chǎn)物CH4,C2H4和C2H2的濃度分布如圖4所示。由圖可看出,反應器內(nèi)溫度由于裂解反應而迅速下降,并在冷卻水的影響下進一步降低,不同實驗條件下所得的淬冷前溫度分別為1 752,1 857 及2 415 K。非均相模型計算所得的淬冷前溫度高于均相模型的模擬值,這是由于焦的生成過程為放熱反應,使得反應段的溫度升高。反應器內(nèi)CH4與C2H4濃度在0.02~0.10mm處達到最大值,隨后便急劇下降,其變化趨勢與均相模型的計算結(jié)果大致相同。C2H2濃度在反應器內(nèi)快速到達最大值后有較為明顯的下降。隨著輸入功率的增加,C2H2濃度達到最大值的速率加快。輸入功率為614.9 kW時,C2H2濃度在112mm(反應的停留時間為1.12 ms)處達到最大,為14.81%;輸入功率為692.8 kW時,在55mm(反應的停留時間為0.52 ms)處達到最大,為14.80%;輸入功率為794.2 kW時,在33mm(反應的停留時間為0.28 ms)處達到最大,為14.43%。由此可見,C2H2最高濃度隨著輸入功率的增加而降低,這是由于反應器內(nèi)的溫度提高加快了固體焦的生成速率,不利于 C2H2濃度的提升。同時可以看到,當停留時間進一步增加,反應更傾向于焦的生成過程,因此減少反應的停留時間能夠有效提升C2H2收率??梢娡ㄟ^縮短反應器長度,降低反應的停留時間至1.0 ms以下能有效提升C2H2收率(約10%左右)。對比中試實驗結(jié)果,如表3所示,非均相模型模擬所得的C2H2產(chǎn)率仍比實際值偏大,但是在實驗過程中,淬冷水與產(chǎn)品氣之間發(fā)生反應所生成的CO甚至CO2會降低C2H2收率??紤]到產(chǎn)品氣中C主要以C2H2的形式存在,可以認為實驗中所檢測到的CO甚至CO2主要由淬冷水與C2H2反應生成。將實驗結(jié)果中C2H2,CO和CO2的收率視為理想淬冷狀態(tài)下的C2H2收率,并考慮氣相產(chǎn)物中的碳收率,非均相模型的模擬值與實驗值有較好的一致性。

        圖4 非均相反應動力學模型模擬結(jié)果中CH4,C2H4與C2H2的濃度分布與溫度分布Fig.4 Distributions of CH4, C2H4and C2H2concentrations and temperature from the simulation of heterogeneous reaction kinetic model

        表3 實驗值與模擬值的比較Table 3 Comparison of experimental data and simulation results

        2.3均相反應動力學模型與非均相反應動力學模型對比

        當輸入功率為692.8 kW時,均相反應動力學模型與非均相反應動力學模型模擬所得的反應器內(nèi)C2H2的濃度分布如圖5所示。由圖可看出,兩種模型給出的反應器內(nèi)C2H2濃度分布在反應段前端0~50mm(反應停留時間為0.5 ms)處幾乎完全相同,隨著停留時間的進一步加長(0.5~4.0 ms),結(jié)焦反應對于裂解氣的組成產(chǎn)生影響。由此可以看出,焦的生成速率遠遠低于氣相反應過程。這意味著當?shù)入x子體反應器的停留時間足夠短時,均相反應動力學模型亦可用于裂解過程的模擬計算。

        圖5 均相模型與非均相模型的乙炔濃度分布Fig.5 Concentration distribution of C2H2from homogeneous and heterogeneous kinetic model under Exp 2

        3 結(jié) 論

        將旋轉(zhuǎn)弧等離子體反應器視為一維平推流反應器模型網(wǎng)絡,結(jié)合裂解反應動力學模型與反應器流動模型,采用CHEMKIN-PRO對丙烷的裂解過程進行數(shù)值模擬。結(jié)果顯示,非均相反應動力學模型的模擬結(jié)果與實驗結(jié)果有更好的一致性。分析反應器內(nèi)產(chǎn)物分布可看出C2H2濃度隨停留時間的增加存在最佳值。當輸入功率為614.9 kW時,C2H2濃度在停留時間為1.12 ms處達到最大,為14.81%;當輸入功率為692.8 kW時,C2H2濃度在停留時間為0.52 ms處達到最大,為14.80%;當輸入功率為794.2 kW時,C2H2濃度在停留時間為0.28 ms處達到最大,為14.43%。隨著停留時間的增加,結(jié)焦過程對裂解氣組成的影響逐漸增大??梢钥闯鼋沟男纬蛇^程在整個裂解過程中是一個動力學限制過程,通過縮短反應器長度,降低反應的停留時間至1.0 ms以下能有效提升C2H2收率(10%左右)。隨著輸入功率增加,反應段內(nèi)C2H2濃度的最大值降低,過高的功率會加速反應體系內(nèi)的結(jié)焦過程,不利于C2H2收率的提升。

        [1]Yan B H, Xu P C, Guo C Y, et al. Experimental study on coal pyrolysis to acetylene in thermal plasma reactors[J]. Chemical Engineering Journal, 2012, 207(SI):109-116.

        [2]Yoon S J, Lee J G. Hydrogen-rich syngas production through coal and charcoal gasification using microwave steam and air plasma torch[J]. International Journal of Hydrogen Energy, 2012, 37(22):17093-17100.

        [3]Hsieh L T, Lee W J, Chen C Y, et al. Converting methane by using an RF plasma reactor[J]. Plasma Chemistry and Plasma Processing,1998, 18(2):215-239.

        [4]Yang Y C, Lee B J, Chun Y N. Characteristics of methane reforming using gliding arc reactor[J]. Energy, 2009, 34(2):172-177.

        [5]Zhao Z L, Huang H T, Wu C Z, et al. Biomass pyrolysis in an argon/hydrogen plasma reactor[J]. Chemical Engineing & Technology,2001, 24(11):197-199.

        [6]Slovetskii D I, Mankelevich Y A, Slovetskii S D, et al. Mathematical modeling of the plasma-chemical pyrolysis of methane[J]. High Energy Chemistry, 2002, 36(1):44-52.

        [7]Fincke J R, Anderson R P, Hyde T, et al. Plasma thermal conversion of methane to acetylene[J]. Plasma Chemistry and Plasma Processing, 2002, 22(1):105-136.

        [8]Tian Y J, Xie K C, Zhu S Y, et al. Simulation of coal pyrolysis in plasma jet by CPD model[J]. Energy & Fuels, 2001, 15(6):1354-1358.

        [9]Yan B H, Xu P C, Jin Y, et al. Understanding coal/hydrocarbons pyrolysis in thermal plasma reactors by thermodynamic analysis[J]. Chemical Engineering Science, 2012, 84:31-39.

        [10]Holmen A, Rokstad O A, Solbakken A. High-temperature pyrolysis of hydrocarbons 1:methane to acetylene[J]. Industrial & Engineering Chemistry Process Design and Development, 1976, 15(3):439-444.

        [11]Fincke J R, Anderson R P, Hyde T A, et al. Plasma pyrolysis of methane to hydrogen and carbon black[J]. Industrial & Engineering Chemistry Research, 2002, 41(6):1425-1435.

        [12]Wang H, Frenklach M. A detailed kinetic modeling study of aromatics formation in laminar premixed acetylene and ethylene flames[J]. Combustion and Flame, 1997, 110(1/2):173-221.

        [13]Frenklach M, Wang H. Detailed mechanism and modeling of soot particle formation[M]. Bockhorn H. Berlin:Springer Berlin Heidelberg,1994:165-192.

        [14]Wen J Z, Thomson M J, Park S H, et al. Study of soot growth in a plug flow reactor using a moving sectional model[J]. Proceedings of the Combustion Institute, 2005, 30(1):1477-1484.

        [15]Wang H, Dames E, Sirjean B, et al. High-temperature chemical kinetic model of n-alkane (up to n-dodecane), cyclohexane, and methyl-,ethyl-, n-propyl and n-butyl-cyclohexane oxidation at high temperatures[EB/OL]. [2010-09-19]. http://web.stanford.edu/group/haiwanglab/ JetSurF/JetSurF2.0/index.html.

        [16]Wang H, You X, Joshi A V, et al. High-temperature combustion reaction model of H2/CO/C1-C4compounds, USC Mech Version II[EB/OL]. [2007-05]. http://ignis.usc.edu/USC_Mech_II.htm.

        [17]Dworkin S B, Zhang Q A, Thomson M J, et al. Application of an enhanced PAH growth model to soot formation in a laminar coflow ethylene/air diffusion flame[J]. Combustion and Flame, 2011, 158(9):1682-1695.

        [18]Sanchez N E, Callejas A, Millera A, et al. Formation of PAH and soot during acetylene pyrolysis at different gas residence times and reaction temperatures[J]. Energy, 2012, 43(1):30-36.

        [19]Indarto A. Soot growth mechanisms from polyynes[J]. Environmental Engineering Science, 2009, 26(12):1685-1691.

        [20]Frenklach M. Reaction mechanism of soot formation in flames[J]. Physical Chemistry Chemical Physics, 2002,4(11):2028-2037.

        [21]Sabbah H, Biennier L, Klippenstein S J, et al. Exploring the role of PAHs in the formation of soot:pyrene dimerization[J]. Journal of Physical Chemistry Letters, 2010, 19(1):2962-2967.

        [22]Wen J Z, Thomson M J, Lightstone M F, et al. Detailed kinetic modeling of carbonaceous nanoparticle inception and surface growth during the pyrolysis of C6H6behind shock waves[J]. Energy & Fuels, 2006, 20(2):547-559.

        Numerical Simulation of Propane Pyrolysis to Acetylene by Thermal Plasma

        Su Baogen, Wu Jianhua, Wen Guangdong, Ma Jie, Fang Jianwei, Xing Huabin, Ren Qilong
        Key Laboratory of Biomass Chemical Engineering of Ministry of Education,Department of Chemical and Biological Engineering,Zhejiang University, Hangzhou 310027, China

        In order to rationally optimize the residence time of propane pyrolysis process in a thermal plasma reactor, the rotating arc plasma reactor was regarded as one-dimensional plug flow reactor network model to analyze the concentration and temperature distributions of the propane pyrolysis in the thermal plasma reactor using the pyrolysis reaction kinetic model and the flow model of the reactor to mathematically simulate the pyrolysis process of propane by CHEMKIN-PRO software. The plasma pyrolysis process of propane was investigated from homogeneous and heterogeneous reaction models, respectively. The simulation results indicated that the heterogeneous reaction kinetic model with the soot formation model was in agreement with the experimental results. With the increase of the reactor length, there was an optimum concentration point of acetylene and the yield of acetylene could be efficiently improved by reducing the residence time of plasma reactor to 1.0 ms.

        numerical simulation; thermal plasma; propane; pyrolysis; acetylene

        TQ018;TQ221

        A

        1001—7631 ( 2016 ) 04—0341—07

        2015-03-19;

        2016-07-13。

        蘇寶根(1975—),男,副教授;任其龍(1959—),男,教授,通訊聯(lián)系人。E-mail:renql@zju.edu.cn。

        猜你喜歡
        丙烷等離子體反應器
        聚丙烯環(huán)管反應器升密操作與控制
        云南化工(2021年11期)2022-01-12 06:06:36
        連續(xù)磁活動對等離子體層演化的影響
        基于低溫等離子體修飾的PET/PVC浮選分離
        流化床丙烷脫氫反應段的模擬及優(yōu)化
        EGSB反應器處理阿維菌素廢水
        等離子體種子處理技術介紹
        上旋流厭氧反應器在造紙廢水處理中的應用
        狀態(tài)監(jiān)測技術在丙烷壓縮機上的應用
        費托合成微反應器研究進展
        化工進展(2015年6期)2015-11-13 00:27:28
        用丙烷作運輸燃料對CNG構(gòu)成了挑戰(zhàn)
        84pao强力打造免费视频34| 国产精品久久久久久久久电影网| av鲁丝一区鲁丝二区鲁丝三区| 国产乱理伦片在线观看| 人妻少妇人人丰满视频网站| 高清中文字幕一区二区三区| 国产综合色在线精品| 伊人久久无码中文字幕| 午夜无码熟熟妇丰满人妻| 蜜臀av在线一区二区尤物| 日韩av无码中文无码电影| 亚洲一区二区三区国产精华液 | 国产性感丝袜在线观看| 天天爽天天爽夜夜爽毛片| 欧美激情二区| 粗一硬一长一进一爽一a视频| 日本一区二区视频免费在线看| 国产98在线 | 日韩| 极品 在线 视频 大陆 国产| 中文字幕亚洲中文第一| 日本精品少妇一区二区三区| 乱码午夜-极国产极内射| 国内精品九九久久精品小草| 狼人伊人影院在线观看国产| 最新系列国产专区|亚洲国产| 国产一级片毛片| 亚洲国产精品嫩草影院久久av| 夜夜躁日日躁狠狠久久av| 精品性高朝久久久久久久| 亚洲精品无人区一区二区三区| 亚洲av综合av一区| 夜先锋av资源网站| 国产成人精品三上悠亚久久| 午夜理论片日本中文在线| 又色又爽又高潮免费视频观看| 国产亚洲精品看片在线观看| 日本精品啪啪一区二区| 午夜免费电影| 黑人巨大videos极度另类| 精品一区二区三区女同免费| 成人免费播放视频777777|