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

催化裂化裝置反再系統動態模擬精細化與控制系統“工藝優先”配對設計

2022-03-03 05:53:06張興碩羅雄麟許鋒
化工學報 2022年2期
關鍵詞:催化裂化模型系統

張興碩,羅雄麟,許鋒

(中國石油大學(北京)信息科學與工程學院自動化系,北京 102249)

引 言

催化裂化在我國的煉油工藝中占有舉足輕重的地位[1],而催化裂化反應-再生系統是煉油流程中重要的二次加工裝置。為了建立自動化優化平臺,需要對煉油裝置進行優化控制,在實際生產中,設計控制回路需要逐一測試然后逐一閉環。為確保安全性和經濟效益,測試控制系統的效果不能用已投入生產的裝置,而未投入使用的裝置由于安全性問題,也不能盲目用于控制系統設計的測試。因此,必須要有一個準確的模型來模擬實際工廠的運行,并能夠進行動態模擬仿真。這種仿真平臺的搭建提供了極大的便利條件,可以在不需要考慮安全問題和經濟損失的情況下進行控制系統設計。

催化裂化反應動力學模型的研究與開發不僅能夠對反應器的設計、產品分布及性質的預測、現場操作的優化、催化劑的選用等提供指導,而且還可加深對反應機理及過程的認識和理解,從而進一步推動工藝技術的發展?;诖呋鸦跓捰凸に囍械闹匾匚唬呋鸦磻獎恿W的研究受到極大重視和關注。隨著反應工程學和計算機技術的發展,出現了集總方法,Weekman[2]最早將集總理論成功運用于催化裂化過程,建立了催化裂化三集總反應動力學模型。之后提出的四集總、五集總等[3-7]都是遵循了三集總反應動力學模型按餾程劃分集總的思想。

在實際化工生產過程中,由于組分的流動需要一定時間,因此信息的傳遞往往存在時滯現象[8],對于大型生產裝置來說這種現象更為明顯,而現有的理論研究大部分都未考慮時滯問題或者簡單地將時滯代入容器的出口、入口處,這就導致模型與實際情況產生偏差。

另外,在氣固相并存的反應器中,氣相的變化與固相不盡相同,由于氣相快速、靈敏的動態特性,氣相總是更快產生變化[9],而由于固相的存在,氣相最終會跟隨固相特性。在大部分學者的研究中都會為了簡化模型,將氣相擬穩態處理,然而這么做的結果就是反應初期的氣相變化與實際不符,且模型中壓力的變化由溫度主導,產生偏差。

本文分別搭建了催化裂化反應器和再生器模型,恢復了催化裂化裝置反應-再生系統的氣相動態和分布參數系統的時滯,進行仿真并分析曲線變化趨勢,驗證了該精細化模型的有效性,并且作為模擬平臺,為后續的控制系統設計提供更接近實際生產情況的催化裂化反再系統模型。

在針對特定對象或過程進行控制系統設計時,過去大部分化工生產過程的控制系統,是根據化工工藝和經驗直接確定的控制結構。在Bristol[10]提出了相對增益陣(RGA)方法后,隨著控制系統理論研究的不斷深入和發展,出現了諸如相對正則化增益陣(RNGA)、相對能量增益陣(REGA)、動態相對增益陣(Dynamic-RGA)等改進方法[11-20]。

與大部分研究中單獨使用化工工藝配對法或系統設計控制配對法不同,本文在進行控制系統設計時,采用“工藝優先”配對,對結構進行多角度分析,先通過工藝配對保證系統穩定,將系統維數降低,然后通過系統理論將剩余變量配對。

1 問題提出

在進行控制系統設計時,大部分研究是單獨采用化工工藝法[21-23]設計或者使用系統理論對模型進行理論設計?;すに嚪梢愿鶕磻^程和以往經驗直接確定控制回路,系統理論對于變量的配對更加嚴謹,并且對相關專業知識的要求較低,應用面更為廣泛。然而對于實際生產過程來說,大多數都是多變量不穩定系統,其中不乏耦合、循環系統,因此,單一的化工工藝確定的控制回路較少,較難分析多變量之間的關聯程度和影響。而單一的系統理論變量配對在對復雜系統進行設計時較為煩瑣,且對大部分實際生產中的不穩定耦合系統的控制回路設計過程也存在諸如傳函求取困難、無法取得穩態特性等問題。

因此本文將結合工藝配對法和系統理論設計配對法,通過工藝優先配對法進行控制系統設計,這種設計方法結合了工藝法的簡便性和系統理論設計的精確性,較為簡便且控制效果良好。這種全局思考將各問題結合分析,在建立精細化模型的同時,在其基礎上進行控制系統的設計。這種綜合協調將需求轉化為實際系統的全部技術工作,對所有系統都具有普遍意義[24]。

本文擬定三級控制回路“工藝優先”設計法。

(1) 以工藝設計確定一級物料平衡控制回路。該控制回路主要保證系統運行以及物料反應過程的穩定性。

(2) 以工藝設計確定二級能量平衡控制回路。該控制回路主要保證系統的長時間運行以及溫度穩定性。

(3) 以系統理論設計確定三級產品質量控制回路。該控制回路影響最小,目的是在一二級控制回路的作用下,對潛在被控變量進行控制,綜合考慮控制回路設計后的復雜度來決定該回路是否有設計需要。

這種“工藝優先”設計法,保證了系統的穩定性,控制了耦合系統的重要變量。這種采用化工工藝法與系統理論設計的相對能量增益矩陣(REGA)相結合的設計方法,在確定變量間關聯和控制系統結構時,可以適配實際生產中多數的不穩定、耦合的復雜系統。

控制系統設計的目的就是對某過程實現控制,然而控制系統不能在實際的生產裝置上進行測試,這時就需要建立過程的數學模型。而現有的文獻對于催化裂化反再系統的研究,大部分都是以氣相穩態且忽略分布參數系統時滯為前提,這樣做導致的結果就是模型精度的降低。而基于這樣的簡化模型設計出的控制系統,最終實際的效果大多和理論上的效果存在一定差距[25]。

一般來說,在化工過程的系統中,傳遞普遍存在時滯現象,在模型化分布參數系統時,普遍會將其中的時滯忽略,導致其不能反映實際的化工生產情況。且在催化裂化反應-再生系統中,為簡化模型,大部分研究都會將氣相按擬穩態處理,但是,在實際工業生產過程中,氣相在初始時刻的變化與固相是有一定區別的,因此氣相擬穩態處理會將氣相的動態特性忽略,此時會導致模型精度進一步降低?;谝陨夏P偷牟粶蚀_,導致最終模型和實際生產的運行結果也會存在一定的差值,所做設想如圖1所示。

圖1 在初始時刻進行輸入階躍后氣相的變化Fig.1 Change of gas phase after input step at initial time

在上述簡化模型的基礎上設計控制系統,本身就受到了一定局限,與實際生產情況的差異,導致了最終使用設計出的控制系統時,實際的控制效果不夠理想甚至不能保證安全運行。

對于氣固共存的催化裂化裝置反再系統,由于與氣相物料相關的壓力和氣相溫度變化速度很快,而固相變量變化較慢,二者不在相同時間尺度上,而一般建模過程中都忽略了氣相的快變化過程,將其擬穩態處理。壓力通過理想氣體狀態方程求得,沒有獨立的動態特性。

本文對催化裂化裝置反再系統進行了建模與仿真。首先,針對催化裂化裝置反再系統,通過實際工業裝置的運行狀態和結構,在已有的再生器和反應器模型的基礎上進行精細化建模。

本文中機理模型由一系列偏微分方程組成,在實際過程中,物料和溫度傳遞需要一定時間,因此,將系統存在的時滯進行恢復。對于反再系統的多時間尺度特性,通常情況下都將與氣相物料相關的動態方程進行擬穩態處理,通過降階達到簡化模型的目的,但也因此無法得到氣相過程的真實特性。此時壓力為主的氣相組分完全依賴于溫度,失去了獨立的動態特性,因而在進行控制等研究時,對溫度控制就相當于對壓力進行控制,這就導致了無法進行嚴格的壓力控制。

在實際化工過程中,壓力變化要比溫度變化快得多,原有的擬穩態模型與實際過程相差較大。如果模型與實際過程相差太大,那么模型的動態過程就不能反映實際過程,模擬仿真結果也將會失去對實際化工過程的指導意義,而在這種模型的基礎上進行的控制系統設計會存在偏差和風險。因此,反再系統的精細化建模尤為重要,本文將對模型進行兩個方面的精細化處理:恢復氣相動態和引入分布參數系統時滯。

2 模型的恢復

基于以往研究,在不加入控制器的情況下,想要保持該化工過程的穩定性,需要控制藏量和壓力恒定。因此在進行模型對比時,設定藏量和壓力為定值,且在研究進行中發現,溫度也是保證反應穩定運行的重要條件,這點將在第3節進行討論。

以本課題組[26-28]的五集總反應動力學模型為基礎,動力學模型表示如下

其中,υ為轉化系數;k為反應速率常數,s-1。

為了確定反應過程中某段容器的組分、溫度等信息,數學模型以微分方程的離散化模型為基礎。在原模型中,提升管、燒焦罐以及密相床的氣體均按擬穩態處理,以提升管為例。

原模型

其中,n為離散化分段數;yAi為提升管第i段原料未轉化率,%;ST為空時,s;kAi為第i段中原料裂化反應速率常數,s-1;GCrg2為再生催化劑循環量,kg·s-1;FO為進料質量流量,kg·s-1;Kh為重芳環吸附常數;CA為碳原子質量分數,%;D1為回煉油與新鮮原料的比值;D2為回煉油漿與新鮮原料的比值;δ1和δ2為物料性質有關常數;pra為反應壓力,Pa;?i為第i段催化劑相對活性。

沉降器壓力模型

其中,praf為沉降器壓力,Pa;Traf為沉降器平均溫度,K;Vraf為沉降器容積,m3;R為理想氣體常數,J·mol-1·K-1;yG|X=1、yN|X=1和yD|X=1分別為氣體、汽油和柴油在出口處產率,%;MWgas、MWN、MWD、MWHLO和MWSlurry分別為氣體、汽油、柴油、回煉油和油漿的平均相對分子質量;FWst、FWpre和FWatmr分別為汽提蒸汽、提升管預熱提升蒸汽和霧化蒸汽流量,kg·s-1;pt為分餾塔頂壓力,Pa;Δpfra為分餾塔塔板液體靜壓,Pa。

在壓力模型中,氣相無動態變化,對壓力變化無影響。

2.1 模型精細化

本文作者在建模時受當時的計算機能力所限,將氣相按擬穩態處理。如今計算機的計算能力足夠支撐各種復雜模型的計算。

針對催化裂化裝置模型的改進,是以前期研究[26-29]為基礎。在原氣相穩態模型的基礎上,不再忽略氣相導數項,恢復氣相動態后,以原料未轉化率為例

原模型中忽略氣相動態導致壓力主要受溫度影響,此時恢復氣相動態后的壓力將受到氣相動態影響。

受以下泛函微分方程模型[30]的啟發,將分布參數系統中物質成分、性質和溫度等傳遞導致的時滯進行還原。

在分布參數系統中,物料和溫度傳遞需要一定時間。在一階分布參數系統中,時滯為兩點之間對流所需時間,而在反應-再生系統的二階分布參數系統中,時滯存在于對流和擴散中,為兩點向中間傳遞所需時間。

在催化裂化反應-再生系統中,擴散所需時間與反應過程相比,為較大的時間尺度。在僅考慮時間尺度時,擴散對系統時滯的取值會有很大影響,然而在同時存在對流與擴散的反應中,擴散時間對反應的影響可以忽略不計,因此在恢復時滯時,忽略在大時間尺度上擴散的影響,在反應過程的短時間內,擴散導致的時滯亦不必考慮。

以提升管一階分布參數系統的原料未轉化率為例,恢復時滯后的模型為

其中,τAi為原料流動的時滯,s。由對流速度及每段提升管長度求出,每段反應的原料未轉化率相對入口處都存在時滯,該時滯為第一段到該段的時滯總和。為突出受時滯影響的變量,其余時變變量在此不進行標識。

根據以上方法,得到恢復氣相動態以及時滯的精細化模型

模型精細化按此類方法逐一進行,在此不再贅述。

2.2 模型對比

保持主風總量不變,保持其他參數不變,對主風量進行階躍測試,即燒焦罐主風增加100 m3·h-1,同時二密相床主風降低100 m3·h-1,得到模型運行曲線如圖2所示。

以汽油產率為例,給出四種對比模型如下。

原模型

分布參數時滯恢復模型

表1 模型的部分相關參數Table 1 Some relevant parameters of the model

其中,yNi和yDi分別為提升管第i段汽油和柴油的產率,%;kNi和kDi分別為第i段中汽油和柴油的裂化反應速率常數,s-1;υAN和υDN分別為原料裂化汽油和柴油裂化汽油的轉化系數;τNi為上一段到當前段反應的汽油對流所需時間,s。

與原模型相比,根據圖2可以看出:

(1)恢復分布參數時滯的模型會在變動產生的約10 s后產生變化;

(2)恢復氣相動態的模型在輸入量變動后,汽油產率和氣體產率等氣相變量,在短時間時滯后開始增加,在約45 s后,與固相動態特性變化趨同;

(3)恢復氣相動態和分布參數時滯的模型,在輸入量變動后,經歷一個短時間時滯,隨后氣相由于自身的快時變特性產生瞬間超出固相的波動變化,在氣相的瞬間波動后,由固相主導的反應特性使得氣相與固相趨同,由于裝置特性,在達到穩定前會一直循環疊加時滯,因此最后會導致大時滯的現象產生,從圖2中可看出,該時滯約為400 s。

圖2 四種模型在控制催化劑藏量和壓力的情況下進行主風量階躍后的曲線Fig.2 The curve for four models with main air volume step under the control of catalyst inventory and pressure

同時,根據燒焦罐一氧化碳氧含量進行輔助驗證,如圖3 所示,燒焦罐主風量增加,精細化模型中溫度短暫下降,由于助燃劑未變化,此時一氧化碳含量會有一個迅速但并不是瞬態的增加,在溫度開始上升后,一氧化碳含量隨之逐漸降低直至穩定。

圖3 燒焦罐主風階躍后一氧化碳含量變化Fig.3 Change of CO content after main air step of regenerator

精細化模型與常用模型最大的區別在于初始短時間的瞬間響應與長時間上的時滯現象,而兩個時滯現象以及氣固相瞬態特性的變化更加接近實際的催化裂化裝置的化工生產過程。

基于以上催化裂化反應-再生系統氣固全動態恢復分布參數時滯的模型,進行控制系統設計。

3 控制系統設計

基于以上模型,進行控制系統設計。

通過分析可以得出對化工生產過程有意義的被控變量和可操縱變量如下。

被控變量:原料未轉化率、柴油產率、汽油產率、氣體產率、燒焦罐催化劑藏量、汽提段催化劑藏量、二密相床催化劑藏量、反應溫度、燒焦罐溫升、沉降器壓力、兩器壓差、再生器壓力、再生煙氣氧含量。

可操縱變量:再生滑閥開度、待生滑閥開度、再生器內循環滑閥開度、富氣壓縮機轉速、再生煙氣雙動滑閥開度、原料預熱溫度、燒焦罐主風量、二密相床主風量、一氧化碳助燃劑投放比。

根據化工過程工藝法,可以直接判斷出一些被控變量之間具有的強關聯性,因此將被控變量簡化:原料未轉化率、柴油產率、汽提段催化劑藏量、二密相床催化劑藏量、反應溫度、燒焦罐溫升、沉降器壓力、兩器壓差、再生煙氣氧含量。

因此得到反再系統傳遞函數為:

其中,G0為9×9 傳遞函數矩陣;Y0為上述簡化后被控變量;U0為上述可操縱變量。

反再系統中,保證化工生產穩定的基礎為催化劑藏量的穩定。在反應過程中,壓力與溫度同時對反應造成影響,但在氣相擬穩態模型中,壓力主要跟隨溫度變化。而恢復氣相動態后,壓力反應變得靈敏,因此壓力控制也是保證反應穩定的重要控制回路。

在第2 節進行實驗時發現,理想控制狀態下控制藏量而不控制溫度時,在約80 min 后,提升管溫度從495℃開始下降,同時燒焦罐溫度從695℃開始下降,會因失溫最終導致反應停止,如圖4 所示,因此在進行控制系統設計時,需要首先對特定變量進行控制以確保反應穩定進行。

圖4 反再系統在不控制溫度時的溫度曲線Fig.4 Temperature curves of reactor-regenerator system without control of temperature

根據化工生產過程的工藝法及以往研究,首先確定四個一級控制回路:

(1) 再生器內循環滑閥-二密相床催化劑藏量控制;

(2)待生滑閥-汽提段催化劑藏量控制;

(3)富氣壓縮機轉速-沉降器壓力控制;

(4)再生煙氣雙動滑閥-兩器壓差控制。

在模型改變后,驗證以上四組變量配對,通過控制實驗發現控制回路依然成立。

通過計算,得到系統的相對能量增益矩陣。由表2可知,剩余變量已經無法繼續進行配對,這也與實際化工生產過程的控制結構相符。因此,得到最終的控制結構變量配對如表3所示。

表2 基于工藝優先的相對能量增益矩陣Table 2 Results of REGA based on prior process analysis

表3 變量配對結果Table 3 Results of variable pairs

確定了以上四個控制回路后,一旦反應向溫度降低的方向偏離,系統會逐漸降溫直至熄火。為保證反應的長久運行,需要保證溫度控制在一定范圍內,因此由化工工藝法確定二級控制回路:再生滑閥-反應溫度控制。

確定了以上控制回路后,需要進行設計的控制系統維數由9×9降低為4×4,大大降低了后續設計的復雜度,之后進行三級控制回路的設計。

將輸入量M序列化,通過平臺搭建的模型得到輸出數據,將輸入輸出數據通過遞推最小二乘法,辨識得到傳遞函數如下

式中,G為4×4 傳遞函數矩陣;Y=[yayd?Trg1yO2]T;U= [Vrg1Vrg2TfreshXs]T。其中,ya為原料未轉化率,%;yd為柴油產率,%;?Trg1為燒焦罐溫升,K;yO2為煙氣氧含量,%;Vrg1為燒焦罐主風量,m3·h-1;Vrg2為二密相主風量,m3·h-1;Tfresh為原料預熱溫度,K;Xs為一氧化碳助燃劑投放比[29,31],%。

由以上控制系統結構分析可知,控制系統設計原則為先控制藏量、氣體和壓力,然后再控制溫度,為了保證化工生產過程“安、穩、優”,需要嚴格按照順序進行。

若由單純的系統理論設計控制結構,有以下幾種方法:

(1) 去除控制回路,求取9×9 多輸入-多輸出系統的傳遞函數然后求得增益矩陣;

(2)去除控制回路,通過階躍輸入獲得穩態輸出與穩態增益然后求得增益矩陣;

(3)在有控制回路的情況下,求取閉環傳函然后得到開環傳函,通過階躍輸入獲得穩態輸出與穩態增益,再取得增益矩陣。

通過該方法求出相對能量增益矩陣,可以看出,最終結果與“工藝優先”配對設計法得到的結果一致,即和傳統控制結構一致。

單純應用以上方法時在耦合、循環的復雜系統中存在一些問題。首先是去掉控制回路后,系統失去了穩定性,而在溫度和壓力的影響下,各輸入-輸出的傳函不夠準確;其次是階躍輸入后在不穩定系統中并不能獲得穩態輸出;另外是在有控制回路的情況下通過開環取得增益矩陣,對于高維多輸入-多輸出的系統來說非常復雜且過程煩瑣,因此在可以使用“工藝優先”設計法的情況下優先考慮“工藝優先”設計法。表4 為直接對高維多輸入-多輸出系統求取的相對能量增益矩陣,其中Tra為提升管反應溫度,K;Praf為沉降器壓力,Pa;dP為兩器壓差,Pa;Wst為汽提段催化劑藏量,t;Wrg2為二密相床催化劑藏量,t;Zcrg2為再生滑閥開度,%;Fsteam為富氣壓縮機氣體流量,t·h-1;Zrgf為再生煙氣雙動滑閥開度,%;Zcst為待生滑閥開度,%;Zcrg21為內循環滑閥開度,%。

表4 相對能量增益矩陣Table 4 Results of REGA

精細化模型與原模型相比有較大變化,為了驗證模型的控制效果,需要與原模型做控制效果對比。雖然控制方法的研究在不斷深入,但由于傳統PID 控制簡單便捷,在實際生產中許多工廠仍在使用[32]。這里使用本課題組所做過研究為參考的原模型PI 控制器參數[33],同時調試出精細化模型的最佳PID 參數。如圖5 所示,以相同工作點為初始值,通過相同的控制回路,以原控制器參數對原模型進行控制仿真,以調試后的控制器參數對精細化模型進行控制仿真,控制器參數如表5 和表6 所示。

圖5 精細化模型的控制效果曲線Fig.5 Control effect curves of refined model

表5 原模型的控制器參數Table 5 Controller parameters of original model

表6 精細化模型的控制器參數Table 6 Controller parameters of refined model

根據圖5可以看出,精細化模型的溫度、壓力和產率曲線有較大波動,而原模型的曲線則相對平穩,顯然原模型的控制效果更好,這也符合了客觀規律,因為原模型在建模時做出了大量假設,精細化模型更加接近實際的化工生產過程。比如精細化模型中,氣相的變化會在短時間內引起較大波動,同時時滯的存在也增加了控制難度,所以其控制效果不會優于原模型。

4 結 論

在以往的催化裂化反應-再生裝置的模型研究中均忽略了氣相動態和信息傳遞時滯的影響,本文在五集總反應動力學模型和前置燒焦罐模型的基礎上,對氣相動態和分布參數系統的時滯進行了恢復,且通過仿真過程驗證了模型的正確性,提高了催化裂化裝置的精細度。而后在該精細化模型的基礎上,對其進行控制系統的設計,通過工藝配對法和系統設計控制配對法相結合,首先通過先驗知識確定了一級控制回路的特定變量配對,然后根據模型恢復過程中遇到的熄火現象,設計了二級控制回路的變量配對以確保反應過程的長久運行,在經過以上變量配對后降低了控制系統的變量維數,再通過REGA方法確定了剩余變量的配對。最終結果也同時驗證了先驗知識的正確性。原模型忽略了氣相動態和時滯,這樣就會導致實際化工生產與其作參考時的控制器設計存在差異,而精細化模型更加接近實際化工反應過程,在該模型基礎上的設計更加具有參考價值。

猜你喜歡
催化裂化模型系統
一半模型
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
提高催化裂化C4和C5/C6餾分價值的新工藝
催化裂化裝置摻渣比改造后的運行優化
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
3D打印中的模型分割與打包
主站蜘蛛池模板: 免费a级毛片视频| 欧美国产在线看| www.狠狠| 天天躁狠狠躁| 五月婷婷综合色| 精品午夜国产福利观看| 亚洲中文在线看视频一区| 99久久精品无码专区免费| 她的性爱视频| 动漫精品啪啪一区二区三区| 欧美日韩北条麻妃一区二区| 最新日本中文字幕| 亚洲熟女中文字幕男人总站| 五月天久久婷婷| 高清欧美性猛交XXXX黑人猛交| 伊人蕉久影院| 久久99蜜桃精品久久久久小说| 亚洲日韩图片专区第1页| 九色国产在线| 国产办公室秘书无码精品| 亚洲福利片无码最新在线播放| 54pao国产成人免费视频| 国产成人91精品| 欧美激情第一区| 97在线公开视频| 性做久久久久久久免费看| 亚洲精品视频免费看| 美女免费黄网站| 久久无码高潮喷水| 国产成人精品优优av| 亚洲AV一二三区无码AV蜜桃| 色婷婷综合激情视频免费看| 一级福利视频| 国产精品区网红主播在线观看| 国产99精品视频| 九九香蕉视频| 香蕉视频在线观看www| 尤物特级无码毛片免费| 中文无码伦av中文字幕| 狠狠v日韩v欧美v| 呦女亚洲一区精品| 亚洲区第一页| 九九免费观看全部免费视频| 老色鬼久久亚洲AV综合| 国产精品9| 熟女成人国产精品视频| 国产精品性| 天天做天天爱夜夜爽毛片毛片| 国产成年女人特黄特色毛片免 | 97精品国产高清久久久久蜜芽| 极品性荡少妇一区二区色欲 | 狠狠色综合网| 国产你懂得| 色婷婷在线影院| 日本黄色不卡视频| 一区二区三区在线不卡免费| 91视频区| 一区二区三区在线不卡免费| 久久久久久久蜜桃| 国产日韩欧美精品区性色| 久久五月天国产自| 久久国产精品影院| 欧美啪啪一区| 国产高潮视频在线观看| 久久精品亚洲热综合一区二区| 国产微拍一区| 国产免费人成视频网| 久久6免费视频| 少妇露出福利视频| 91久久精品日日躁夜夜躁欧美| 无码专区国产精品第一页| 亚洲欧洲免费视频| 国产精品99久久久久久董美香| 亚洲欧美不卡中文字幕| 色婷婷狠狠干| 真实国产精品vr专区| 色婷婷狠狠干| 亚洲视频三级| 亚洲视屏在线观看| 色亚洲激情综合精品无码视频| 五月婷婷伊人网| 亚洲成在线观看|