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

委內瑞拉減壓渣油的熱解特性及其動力學研究

2022-09-13 13:05:54李家州郭子千張玉明熊青安
石油學報(石油加工) 2022年5期
關鍵詞:模型

李家州, 喬 沛, 郭子千, 張玉明, 張 煒, 熊青安

(中國石油大學(北京) 重質油國家重點實驗室,北京 102249)

進入新世紀,全球石油資源需求量日益增長,而世界原油重質化和劣質化趨勢不斷增加,供需之間的矛盾不斷加深。全球重油資源量巨大,約占剩余石油資源的70%以上[1]。中國原油品質大多中等偏重,500 ℃以上減壓渣油的產率占40%~50%[2],并且中國石油資源消費巨大,2020年中國原油對外依存度已達到約72%。根據《BP世界能源統計年鑒2020》報道,2019年全球石油消費量日均增加90×104bbl,中國石油消費量日均增加68×104bbl(1 bbl=159 L)[3]。原油資源變重,原油加工過程中渣油的產量相應增加[4]。隨著市場對輕質油品的需求量日益增加以及環保要求日益嚴格,如何實現重油的清潔高效利用已成為當前煉油工業面臨的重大問題。

劣質重油(包括常壓重油、稠油、超稠油以及原油加工過程中產生的減壓渣油和油漿等)通常具有高殘炭值、高碳/氫比、高黏度和密度、易縮合生焦等特性,這些特性對現有的重油加工技術提出了巨大挑戰。重油改質工藝主要分為脫碳和加氫,脫碳工藝主要包括催化裂化、減黏裂化、焦化(延遲焦化、靈活焦化、流化焦化);而加氫工藝則包括加氫裂化、加氫精制、加氫脫硫等。焦化技術由于具有原料適應性強、靈活性好、脫碳率高和投資低等優勢,已成為重油預處理技術和研究開發的熱點[5]。熱解是劣質重油熱化學轉化過程中的重要環節,研究重油熱轉化動力學對于理解熱解過程、確定脫碳工藝參數具有重要意義,通過熱解動力學可以深入了解反應過程及反應機理,預測反應速率、復雜程度和反應難易程度等,可為劣質重油相關工藝的開發和高附加值利用提供依據[6-7]。

常用的動力學分析方法可以分為模型法和無模型法[8]。模型法主要有Coats-Redfern[9]和Dolye等方法,需要先構建動力學模型,若選取不當會造成較大誤差。無模型法(又稱為等轉化率法)主要有Friedman-Reich-Levi(Friedman)[10]、(FWO)[11]、Kissinger-Akahira-Sunose(KAS)[12]等,該方法避免了模型函數的選取,可獲得較為準確的動力學參數。許多研究者傾向于利用一級動力學模型描述劣質重油的熱轉化反應[13-14]。高岱巍等[15]采用熱重分析儀考察華北油漿的熱轉化過程,在溫度為300~596 ℃內采用一級反應動力學方程擬合,計算得到的表觀活化能約為108 kJ/mol。諶家豪等[16]采用一級反應動力學擬合3種催化油漿凈化后的渣漿,其活化能分別為30.40、39.90和32.30 kJ/mol。Che等[17]研究伊朗減壓渣油的熱解行為及動力學,采用FWO法計算得到的活化能為160.00~217.00 kJ/mol。Zhang等[18]將脫油瀝青熱解過程分為兩階段,分別采用一階積分模型計算活化能和指前因子,第一階段225~425 ℃范圍內活化能為37.69~47.91 kJ/mol,第二階段425~528 ℃范圍內活化能為110.00~121.33 kJ/mol。廖澤文等[19]對瀝青質的熱解動力學進行研究,采用Friedman法計算得到的活化能為159.00~455.00 kJ/mol。等轉化率法將劣質重油熱解反應視作單一組分反應,對于精確的反應模型和動力學參數分析可能會有一些系統誤差。劣質重油是多種化合物的混合物,其熱解反應是由多種反應綜合構成的復雜過程,因此將其視為單一組分反應不能夠準確描述其熱解特性[20-22]。分布活化能模型(DAEM)法作為一種高精度的分析方法,其假設有許多具有不同動力學參數的若干一階或n階不可逆平行反應同時發生,通過對其若干偽組分動力學參數進行計算,如平均活化能、指前因子等,實現對復雜反應的分析[23]。DAEM法目前已經被應用于生物質[24-25]、煤炭[26-27]、油砂瀝青[28-29]等多種復雜熱解反應,但是對于減壓渣油的相關研究還在探索過程中。高斯型DAEM法動力學參數之間有著很強的補償效應,通常為了獲得較好的擬合效果,需要對指前因子(A)采用試錯法進行選擇,從而導致表觀活化能(E0)呈現出非唯一解的現象,即所謂的動力學補償效應,因此通常將指前因子A視作常數,從而保證動力學參數的唯一解[30-31]。

筆者采用熱重-質譜(TG-MS)聯用技術研究減壓渣油的熱解反應特性以及氣體產物釋放規律。TG-MS聯用技術已在生物質和煤熱轉化研究和利用中得到廣泛的應用[32-38],而在減壓渣油熱轉化中的應用則鮮有報道。基于熱重數據,利用Friedman、FWO和KAS 3種等轉化率方法對減壓渣油的熱解反應過程進行動力學分析,求取其熱解動力學參數,并比較3種等轉化率方法在減壓渣油熱解過程中的適用性;基于3種無模型動力學方法的結果,進一步使用分布活化能模型DAEM法對升溫速率10 ℃/min條件下的熱重分析(TG-DTG)數據進行擬合,獲得與實驗數據擬合程度較高的四組分高斯曲線及動力學參數,并通過計算獲得活化能分布曲線,推測得到總活化能分布范圍。

1 實驗部分

1.1 實驗原料

選取委內瑞拉減壓渣油作為實驗原料,由中國石油石油化工研究院提供,原料性質如表1所示。由表1可知:減壓渣油在100 ℃時的運動黏度為19647.23 mm2/s,流動性很差;其康氏殘炭值為21.15%,氫/碳原子比較低(1.41),這說明減壓渣油的芳香結構含量較高,芳香環系的稠合程度較高,生焦傾向較大。膠質和瀝青質含量所占比例較高,分別達到質量分數24.70%和12.41%,兩者是渣油熱裂解生焦的主要來源。

表1 委內瑞拉減壓渣油的基本性質Table 1 Basic properties of Venezuela vacuum residue

1.2 實驗裝置及方法

采用日本SⅡ TG/DTA 7300型熱重分析儀考察委內瑞拉減壓渣油的熱解特性。實驗以N2作為惰性保護氣和載氣,N2純度≥99.999%,其體積流量為80 mL/min,設定溫度從30 ℃升高至800 ℃,采用的升溫速率分別為5、10、15和20 ℃/min。實驗中每次稱取樣品質量為4 mg左右,熱天平自動記錄樣品質量的變化信號。實驗開始前進行3次重復實驗,保證數據相對誤差小于5%,選取3次實驗數據的平均值計算動力學參數。同時將熱重分析儀和質譜儀(德國ThermoStar GSD 320型)聯用,獲取實時的氣體釋放信息。

1.3 減壓渣油熱解反應動力學分析

委內瑞拉減壓渣油熱解過程中質量損失反應轉化率按照式(1)計算。

(1)

式(1)中:x為轉化率;m0、mt和mf分別為樣品初始質量、t時刻質量和熱解終止質量,mg。

委內瑞拉減壓渣油熱解過程中反應速率方程按照式(2)計算。

(2)

式(2)中:t為委內瑞拉減壓渣油熱解時間,min;β為升溫速率,℃/min;T為渣油熱解溫度,℃;k為反應速率常數,min-1;f(x)為反應機理函數。

反應速率常數k遵循Arrhenius定律,為:

(3)

式(3)中:A為指前因子,min-1;E為表觀活化能,kJ/mol;R為氣體常數,8.314×10-3kJ/(mol·℃)。

通常,動力學機理函數f(x)的積分形式為:

(4)

由式(2)~式(4)可得:

(5)

通過式(2)~式(5)可得以下通用公式[39-40]:

Friedman法:

(6)

FWO法:

(7)

KAS法:

(8)

根據上述公式,對熱重數據進行處理,分別以ln(dx/dt)、lnβ和ln(β/T2)對1/T進行線性擬合,根據斜率計算活化能。筆者將委內瑞拉減壓渣油的熱解反應機理函數視為一級反應[41],即f(x)=1-x,相應的G(x)=-ln(1-x),再結合擬合曲線截距,可求得指前因子。

一階反應DAEM模型的方程式[24,42]如式(9)所示。

(9)

其中f(E)為反應活化能分布函數。

此時DAEM模型的待求變量實質上是A和f(E),其中常用的f(E)函數有高斯分布、logistic分布、指數分布和Weibull分布,采用最常用的高斯分布函數,分布函數如式(10)所示。

(10)

式(10)中:σ為正態分布里的標準差;E0為反應的平均活化能,kJ/mol。

為了減少動力學補償效應的作用,其中指前因子A通常確定為固定值,因此可以得到DAEM的導數方程如式(11)所示。

(11)

式(11)中:T0為初始時刻的絕對溫度,℃。

采用模式搜索法進行擬合的過程中,隨著迭代次數的增加,當殘差值(Residual error,Re)(見式(12))不再明顯變化時即認為達到最優值,設定Fit作為評價模型擬合程度的擬合度偏差值(見式(13)),該值越低說明動力學參數的擬合效果越好。

(12)

(13)

式(12)和式(13)中:i為研究的數據點;n為數據點的個數;(dx/dT)max表示實驗求出的轉化率導數的最大值。

2 結果與討論

2.1 委內瑞拉減壓渣油熱重分析

2.1.1 熱重曲線分析

委內瑞拉減壓渣油在升溫速率為10 ℃/min下的TG和DTG曲線如圖1所示。從圖1可知,委內瑞拉減壓渣油的熱解反應溫度區間較寬。由于減壓渣油的沸點較高,因此在低于179 ℃范圍較少發生熱解,主要為吸附水分和氣體的脫附;在179~300 ℃之間的質量損失較小,主要為渣油中的輕質組分裂解為低分子產物;在300~490 ℃之間的質量損失較大,此階段包含基于自由基反應機理的減壓渣油熱轉化反應,主要包括環烷烴和芳烴的側鏈斷裂,環烷烴的脫氫和開環,芳香烴的脫氫縮合,烴類中C—C鍵、C—H鍵和非烴類中C—雜原子鍵(如S、N)的斷裂等[43]。高于490 ℃范圍的TG和DTG曲線變化較為平緩,此階段為炭化階段,稠環芳烴脫氫縮合成焦炭,繼續升溫只有少量氣體逸出導致質量損失的發生。減壓渣油在升溫速率10 ℃/min時熱解主要反應溫度段為179~490 ℃,總質量損失率為77.54%,DTG中的質量損失峰在446 ℃達到最大,最大質量損失速率為317.38 μg/min。

圖1 委內瑞拉減壓渣油的TG和DTG曲線Fig.1 TG and DTG curves of Venezuela vacuum residueHeating rate (β) 10 ℃/min

2.1.2 氣體釋放規律分析

當升溫速率為10 ℃/min時,委內瑞拉減壓渣油的熱解氣體(CH4、H2O、CO和CO2)釋放規律如圖2所示。圖2(a)為CH4的逸出質譜曲線,從400 ℃起CH4析出量迅速增加,并在462 ℃時達到峰值,在500~800 ℃區間CH4析出量逐步降低。減壓渣油中烷烴的斷鏈和芳香烴、環烷烴及非烴類的側鏈斷裂會生成大分子的烷基自由基,部分烷基自由基與脫氫產生的氫自由基反應會生成大量的CH4[44]。圖2(b)為H2O的逸出質譜曲線,從室溫到300 ℃溫度區間逸出的H2O主要為減壓渣油物理吸附的H2O以及C—OH鍵斷裂生成的H2O,其達到峰值的溫度為177 ℃。在400~700 ℃溫度區間內,出現H2O的第2個析出峰,這是因為C—O鍵斷裂,其達到峰值的溫度為463 ℃。由于C—O鍵的鍵能大于C—OH鍵的鍵能,故需要更高的溫度才能斷裂。減壓渣油C—O鍵斷裂生成H2O的強度明顯大于C—OH鍵斷裂生成H2O的強度,故其主要的含氧基團為C—O鍵[45]。圖2(c)為CO的逸出質譜曲線,從熱解初始階段就有大量CO析出。在300~700 ℃溫度區間內由于自由基斷裂和縮合反應劇烈,出現CO的最大析出峰,其達到峰值的溫度為459 ℃,接近DTG曲線的最大質量損失速率溫度。圖2(d)為CO2的逸出質譜曲線,從400 ℃時CO2析出量逐漸增加,其達到峰值的溫度為456 ℃,減壓渣油中羧基熱穩定性差,高溫條件下易發生脫羧反應而生成CO2[46]。

根據上述熱解氣體的實時釋放規律,大部分熱解氣初始產生溫度在400 ℃左右,結合委內瑞拉減壓渣油升溫速率10 ℃/min時的熱重曲線,委內瑞拉減壓渣油在400 ℃時的質量損失率為48.13%,其整體熱解過程的質量損失主要是由于裂化和縮合作用。

圖2 委內瑞拉減壓渣油的氣體逸出質譜曲線Fig.2 Gas escape mass spectrometry curve of Venezuela vacuum residue(a) CH4; (b) H2O; (c) CO; (d) CO2

2.1.3 升溫速率對熱重分析結果的影響

委內瑞拉減壓渣油在升溫速率分別為5、10、15和20 ℃/min下的TG、DTG與溫度之間的關系如圖3所示。從圖3可以看出,升溫速率不同,某一時刻的質量損失量有所不同。由于升溫速率不同,熱量至外向內傳遞的速度就不同,升溫速率較低時,樣品受熱充分,質量損失量大,重油熱解的起始溫度和終止溫度較低[47]。升溫速率越高,達到相同質量損失所對應的熱解溫度越高,DTG峰向高溫方向移動。

圖3 不同升溫速率(β)下委內瑞拉減壓渣油的TG和DTG曲線Fig.3 TG and DTG curves of Venezuela vacuum residue at different heating rates (β)(a) TG; (b) DTG

不同升溫速率下對應的最大質量損失速率(Pmax)及其所對應的溫度(Tmax)如表2所示。從表2可知,隨著升溫速率的增加,最大熱解質量損失速率所對應的溫度升高,產生熱滯后現象[48]。另外,升溫速率高,樣品在特定溫度下停留時間較短,完成整個熱解過程所需要的時間減少。

表2 不同升溫速率(β)下最大質量損失速率(Pmax)及其對應的溫度(Tmax)Table 2 Maximum mass loss rate (Pmax) and its correspondingtemperature (Tmax) at different heating rates (β)

2.1.4 轉化率與溫度的關系

圖4為委內瑞拉減壓渣油在不同升溫速率下熱解轉化率與溫度的關系。從圖4可知,相同轉化率下,反應溫度隨升溫速率的增加而增高。在相同溫度下,升溫速率越小,反應轉化率越大。

圖4 委內瑞拉減壓渣油在不同升溫速率(β)下熱解轉化率(x)與溫度(T)的關系Fig.4 Relationship between pyrolysis conversion rate (x) and temperature (T) of Venezuela vacuum residue at different heating rates (β)

2.2 委內瑞拉減壓渣油熱解反應動力學分析

為了獲得Friedman、FWO和KAS 3種模型的動力學關鍵參數,根據圖4的實驗數據,得到轉化率x分別為0.1、0.2、…、0.8時委內瑞拉減壓渣油對應的溫度(T)、ln(dx/dt)、lnβ和ln(β/T2)值,然后分別以ln(dx/dt)、lnβ和ln(β/T2)對1/T進行擬合,得到不同轉化率下的活化能和指前因子,擬合結果如圖5所示。KAS法和FWO法對式(5)在0到T溫度區間內進行了近似溫度積分,這導致隨轉化率變化的活化能在積分區間內被“平均化”,會產生一定的計算誤差。表3列出了通過Friedman法、FWO法和KAS法求得的動力學參數。由于轉化率高于0.8時的相關系數較低,因此不作考慮[49-50]。從表3可知,在轉化率0.1~0.8范圍內,3種模型線性擬合程度均較高,R2都在0.96以上,計算得到的委內瑞拉減壓渣油平均活化能分別為131.37(Friedman法)、115.25(FWO法)和110.69 kJ/mol(KAS法)。

委內瑞拉減壓渣油通過3種等轉化率方法獲得的活化能與轉化率之間的關系曲線如圖6所示。從圖6可知,減壓渣油的活化能不是恒定值,這是由于重油熱解是一個復雜的多步反應過程,不同的反應需要的能量不同[51]。在熱解過程中,重油四組分活化能由小到大順序為飽和分、芳香分、膠質、瀝青質[52]。在反應初期,劣質重油分子中不穩定的鍵(以富含直鏈烴類的飽和分為主)由于鍵能較低而首先發生斷裂,反應的活化能較小。隨著溫度的升高,熱解程度加深,分子中較強的鍵開始斷裂,生成大量氣體產物,反應的表觀活化能增大。在相同轉化率下,FWO和KAS法計算得到的活化能值相差較小,Friedman法計算得到的活化能值稍微偏大。這是由于FWO和KAS法是一種積分等轉化率法,而Friedman法是一種微分等轉化率法,而數值微分計算會放大實驗噪聲,導致Friedman法對實驗噪聲較為敏感[53]。

圖5 Friedman、FWO和KAS法得到的委內瑞拉減壓渣油擬合曲線Fig.5 Fitting curves of Venezuela vacuum residue obtained by Friedman, FWO and KAS methods(a) Friedman; (b) FWO; (c) KAS

表3 由Friedman、FWO和KAS法得到的活化能(E)、指前因子(A)和相關系數(R2)Table 3 Activation energy (E), pre-exponential factor (A) and correlation coefficient (R2) obtained by Friedman, FWO and KAS methods

委內瑞拉減壓渣油通過3種等轉化率方法獲得的動力學參數的范圍見表4。不同方法計算得到的指前因子A和活化能E之間存在一定差距,這是由于動力學補償效應造成的[54]。為了盡量減少動力學補償效應帶來的誤差,應選擇指前因子和活化能相差不大的方法。由于FWO和KAS法計算得到的動力學參數相差較小,但FWO法在溫度積分時的精確度更高[39],故可采用FWO法求得動力學參數。采用FWO法計算得到減壓渣油在轉化率為0.1~0.8范圍內的活化能為56.77~178.91 kJ/mol。

圖6 3種等轉化率方法獲得委內瑞拉減壓渣油的活化能(E)與轉化率(x)之間的關系Fig.6 Relationship between the activation energy (E) and conversion rate (x) of Venezuela vacuum residue obtained by three iso-conversional methods

減壓渣油主要由飽和分、芳香分、膠質和瀝青質四組分構成,因此假設減壓渣油包含4個未知化學成分,首先采用高斯DAEM模型對減壓渣油在升溫速率10 ℃/min下的TG-DTG數據進行分峰擬合,4個假定組分分別標注為Gauss-1、Gauss-2、Gauss-3、Gauss-4,分別對應飽和分、芳香分、膠質和瀝青質,疊加峰標注為Guass,實驗測量值標注為Experiment,各組分分峰結果如圖7所示。將實驗測量值曲線與4個假定組分疊加曲線進行比較,可以看出二者的擬合程度是非常接近的。

為了減少動力學補償效應的影響,在擬合過程中預先設定飽和分、芳香分、膠質和瀝青質的指前因子分別在1012、109、1015、1024的數量級內,分別對4個假定組分進行擬合。圖8為減壓渣油四組分DAEM模型動力學參數的計算結果。根據圖8中動力學參數可以對4個假定組分各自的計算值曲線(分別標注為Calculation-1、Calculation-2、Calculation-3、Calculation-4)與實驗值曲線(分別標注為Experiment-1、Experiment-2、Experiment-3、Experiment-4)作對比,可以看出各假定組分的計算值與實驗值的重合度較高,擬合效果較好。

表4 3種等轉化率方法獲得的動力學參數Table 4 Kinetic parameters obtained by three iso-conversional methods

圖7 委內瑞拉減壓渣油四組分DAEM模型擬合曲線Fig.7 Four-component DAEM fitting curve of Venezuela vacuum residueHeating rate (β) 10 ℃/min

表5為四組分DAEM模型動力學參數及評價指標。從DAEM擬合結果可以看出,Fit值均小于2,能夠基本滿足擬合要求,且相關系數R2均在0.99以上,這表明DAEM模型對于減壓渣油復雜的熱解過程有著很高的適用程度。其次對計算得到的各假定組分的平均活化能進行對比后發現,4種假定組分呈現出E4>E3>E1>E2的規律,同時計算得到的各組分峰占比也基本符合委內瑞拉減壓渣油4種組分的比例,這也說明了對于采用四組分DAEM模型進行擬合的正確性。

如圖9所示,f(E)-1,f(E)-2,f(E)-3,f(E)-4分別對應四種組分的活化能分布曲線,將其疊加后得到委內瑞拉減壓渣油整體的f(E)。委內瑞拉減壓渣油的活化能主要分布在100~250 kJ/mol之間,并且在100~200 kJ/mol之間的活化能分布概率更高,對4種假定組分各自的活化能進行加權平均求和后,計算得到的委內瑞拉減壓渣油的平均活化能為190.11 kJ/mol。相比于采用Friedman、FWO和KAS 3種等轉化率法計算轉化率為0.1~0.8范圍內的活化能,DAEM法可計算重油(轉化率0~1范圍)的總熱解反應活化能,因此DAEM求得的反應活化能更接近減壓渣油實際的活化能范圍。

圖8 委內瑞拉減壓渣油4個假定組分各自的DAEM模型擬合結果與實驗值的對比Fig.8 Comparison of DAEM fitting results of 4 hypothetical components of Venezuela vacuum residue with experimental values(a) Saturates; (b) Aromatics; (c)Resins; (d) Asphaltenes

f(E)-1, f(E)-2, f(E)-3, f(E)-4 correspond to the activation energydistribution curves of the four components, respectively. Aftersuperimposing obtain the overall Venezuela vacuum residue oil f(E).圖9 委內瑞拉減壓渣油的四組分活化能(E)分布Fig.9 Distribution of four-component activation energy (E) of Venezuela vacuum residue

3 結 論

選取委內瑞拉減壓渣油作為典型的劣質重油原料,采用TG-MS聯用技術研究了不同升溫速率下重油熱解質量損失變化及氣體釋放規律;并采用Friedman、FWO和KAS 3種等轉化率法和DAEM模型法對反應過程進行了動力學計算,求解了表觀活化能及指前因子等動力學參數,所得主要結論如下:

(1)委內瑞拉減壓渣油的熱解反應階段為179~490 ℃,總質量損失率為77.54%,質量損失峰在446 ℃達到最大,最大質量損失速率為317.38 μg/min。

(2)相比Friedman和KAS模型,FWO模型能更好地描述委內瑞拉減壓渣油熱解過程,采用FWO法計算的委內瑞拉減壓渣油在轉化率為0.1~0.8范圍內的熱解活化能為56.77~178.91 kJ/mol。

(3)利用DAEM模型法將委內瑞拉減壓渣油分為4個假定組分進行分峰擬合,同時求取四組分動力學參數,并據此獲得委內瑞拉減壓渣油總活化能分布曲線,在升溫速率10 ℃/min下活化能主要集中在100~250 kJ/mol范圍內,通過加權求和獲得平均活化能為190.11 kJ/mol,這與渣油實際熱反應活化能更加接近。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美一级99在线观看国产| 一本一道波多野结衣av黑人在线| 午夜视频免费一区二区在线看| 99青青青精品视频在线| 亚洲欧美综合另类图片小说区| 看国产一级毛片| 日韩国产一区二区三区无码| 午夜视频在线观看免费网站| 欧美另类精品一区二区三区 | 欧美不卡视频在线观看| 亚洲精品天堂自在久久77| 欧美成人免费一区在线播放| 国产精品13页| 亚洲成肉网| 国产精品xxx| 色AV色 综合网站| 中文无码精品A∨在线观看不卡 | 91麻豆精品视频| 三上悠亚一区二区| 国产超薄肉色丝袜网站| 国产成人免费手机在线观看视频| 波多野结衣在线se| 色哟哟国产精品| 国产在线视频导航| 久久狠狠色噜噜狠狠狠狠97视色| 欧美成a人片在线观看| 五月婷婷激情四射| 国模极品一区二区三区| 国产丰满成熟女性性满足视频| 欧美在线伊人| 国产精品久久久久久久久kt| 亚洲精品视频网| swag国产精品| 国产成人乱无码视频| 国产黄在线免费观看| 欧美国产日韩一区二区三区精品影视| 国产SUV精品一区二区6| 久久综合伊人77777| 亚洲自拍另类| 国内精品自在自线视频香蕉| 幺女国产一级毛片| 国产免费精彩视频| 国产精品视频3p| 亚洲无码日韩一区| 激情综合网激情综合| 亚洲国产成人精品无码区性色| a亚洲天堂| 国产成人亚洲精品无码电影| 亚洲一本大道在线| 国产欧美精品午夜在线播放| 免费AV在线播放观看18禁强制| 91成人在线免费观看| 国产噜噜噜视频在线观看| 青青青草国产| 亚洲swag精品自拍一区| 久久激情影院| 国产精品女主播| 在线欧美日韩| 激情亚洲天堂| 99视频在线免费观看| 亚洲成av人无码综合在线观看| 亚洲成人播放| 久草网视频在线| 久久黄色小视频| 免费在线观看av| 又粗又硬又大又爽免费视频播放| 美女视频黄又黄又免费高清| 欧美一级高清片欧美国产欧美| 亚洲精品国产综合99| 精品成人一区二区三区电影| 欧美一级高清免费a| 有专无码视频| 1769国产精品免费视频| 亚洲无码日韩一区| 亚洲AV人人澡人人双人| 成人91在线| 国产激爽大片在线播放| AV不卡在线永久免费观看| 国产在线精品人成导航| 国产91av在线| 国产在线97| 国产在线精彩视频二区|