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

精細搜索策略應用于質量交換網絡綜合

2022-08-10 09:49:34楊嶺崔國民周志強肖媛
化工學報 2022年7期
關鍵詞:優化結構質量

楊嶺,崔國民,周志強,肖媛

(上海理工大學能源與動力工程學院,上海 市動力工程多相流動與傳熱重點實驗室,上海 200093)

引 言

隨著能源與環保問題日益嚴峻,我國提出了“碳達峰·碳中和”的戰略目標[1],這對企業降低排放的能力提出新的要求,使得發展新的節能減排和清潔生產技術成為當下的研究熱點。質量交換網絡(mass exchange network,MEN)[2]作為一項廣泛應用于化工原料的提純、萃取、吸附和污染物去除技術,可顯著降低廢水、廢氣的產生和質量分離劑(mass separating agent,MSA)的消耗,近年來受到越來越多的關注。

質量交換網絡的概念由El-Halwagi等[3]于1989年提出,是在換熱網絡研究的基礎上發展而來的,因此許多學者[4-6]將換熱網絡的夾點法與數學規劃法運用到質量交換網絡上。El-Halwagi 等[3]率先使用熱力學夾點法設計質量交換網絡,在濃度區間表中確定夾點,以此確定具有最大質量交換量的初始網絡,并根據初始網絡改進為成本最低的網絡。Hallale 等[7-8]提出了一種新的y-y*復合曲線圖,用于處理涉及夾點同一側存在多個MSA 的問題,并指出質量交換單元數量最小并不一定使成本最低,網絡設計時應充分考慮傳質驅動力。Gadalla[9]將夾點分析原理與質量交換網絡聯系起來,繪制富流股濃度與等效貧液流濃度的對比圖,所得圖用于分析當前網絡的性能,并設計新的網絡結構。Yanwarizal等[10]通過構建濃度與質量負荷圖同步確定最小MSA 用量和貧富流股間的傳質負荷。

此外,數學規劃法在質量交換網絡綜合上也有著廣泛的應用。El-Halwagi 等[11]使用兩階段合成方法,第一階段使用線性規劃確定外部MSA 的最低成本和夾點,第二階段使用混合整數線性規劃減少質量交換器的數量,從而獲得綜合費用最低的網絡結構。Isafiade 等[12]提出了一種基于區間的混合整數非線性規劃(mixed integer and non-linear programming,MINLP)超結構模型,可同時利用超結構與夾點技術的優點,此外,還規定了分流必須以等濃度混合,雖然簡化了模型的復雜度,但一定程度上限制了模型的求解域。李紹軍等[13]建立了基于分級超結構的無分流質量交換網絡MINLP 模型,并采用Alopex 算法與遺傳算法結合的混合算法進行求解,可同時權衡操作費用和投資費用。謝會等[14]建立了一種基于區間的無分流超結構模型,其中過程貧流股與外部貧流股具有同等匹配機會,且允許重復匹配產生,并將列隊競爭算法應用于該模型。都健等[15]以傳質濃度差作為優化變量,并考慮了多組分系統,利用改進的遺傳算法進行求解,取得了較好的優化結果。侯創等[16]在數學模型中引入了取整函數,并通過LINGO 求解器求解,解決了理論塔板數和實際塔板數存在差異的問題,提高了塔板的利用效率的同時降低了年綜合費用。

啟發式算法與MINLP 問題的適應性和強大的搜索能力在許多研究中得以驗證,但其難以兼顧全局搜索和局部搜索,一定程度上限制了啟發式算法的求解能力。鑒于此,提出一種新的應用于質量交換網絡的精細搜索策略,分別利用了啟發式方法簡單有效、求解域大的特點和確定性方法精度高、收斂快的特點,建立兩種精細搜索方式,分別為帶有個體回代與分化的高精度強制進化隨機游走算法

(random walk algorithm with compulsive evolution,RWCE),以及坐標輪換法與黃金分割法結合的方法。最后通過兩個算例驗證了兩種精細搜索方法的有效性。

1 數學模型

1.1 問題描述

根據目標組分的濃度不同,將過程流股分為貧流股和富流股。為了使富流股中的組分濃度滿足出口濃度要求,通過質量分離器將富流股中所含目標組分通過吸收、解吸、吸附、萃取、離子交換等方式轉移到貧流股中。質量交換網絡優化的目的是在滿足出口濃度要求以及傳質可行性的前提下,獲得最具經濟性的網絡結構。

圖1 質量交換網絡示意圖Fig.1 Sketch map of mass exchange network

貧富流股之間的傳質存在相平衡關系,如式(1)所示,對于稀溶液,在所涉及濃度范圍內,可假設相平衡關系是線性的,且與富流股中其他可溶性物質無關。

圖2 貧富流股的相平衡關系Fig.2 Phase equilibrium of rich and lean stream

1.2 目標函數

目標函數是獲得質量交換網絡的最低年度綜合費用(total annual cost,TAC),包括質量交換器的投資費用C1和外部貧流股的操作費用C2,如式(2)所示。操作費用取決于外部貧流股的流量,而對于不同工況的分離方式,投資費用計算方式不盡相同。當使用板式塔時,通過增加塔板數來增加質量交換器的傳質能力,此時的投資費用由塔板數決定,而使用填料塔時,投資費用由質量交換器的質量決定。操作費用如式(3)所示,投資費用具體計算方式在算例驗證部分詳細介紹。

1.3 有分流節點非結構模型

以往的研究中常采用分級超結構(stage-wise superstructure, SWS)模型[17]合成HEN、MEN,而本文基于換熱網絡優化的思想,將有分流節點非結構模型(nodes-based non-structural model,NNM)[18]應用于質量交換網絡的優化,該模型以節點表示傳質單元潛在的位置,具有更強的靈活性和更寬廣的求解域。模型如圖3 所示,結構中有NR股富流股和NL股貧流股,箭頭表示流股流動方向,所有質量交換單元均為逆流布置。富/貧流股上分別設置NRG和NLG個分流組,富/貧流股經過每個分流組時,分裂為NRS和NLS個分流,該組所有的分流在末端進行非等濃度混合后流入下一分流組。每個分流分別包含NRD和NLD個節點,因此所有的富流股節點數為NR×NRG×NRS×NRD,貧流股總節點數則為NL×NLG×NLS×NLD。初始狀態下,整個結構為空結構,優化開始后隨機在富流股和貧流股上抽取兩個節點,連接形成質量交換單元。

圖3 有分流節點非結構模型Fig.3 Schematic of the nodes-based non-structural model with streams split

1.4 約束條件

2 強制進化隨機游走算法

本文基于節點非結構模型,采用RWCE 算法[19]同步優化質量交換網絡,該算法已在換熱網絡問題上取得了許多優秀的成果[20-21],驗證了該算法強大的全局搜索能力。RWCE 算法的基本原理為:種群中的個體獨立進化,每個個體從空結構開始,以目標函數減小為進化方向,通過隨機生成傳質單元和隨機增減傳質量、分流比及流量實現整型變量和連續變量的同步進化,同時,有一定的概率接受目標函數較高的結構,具備跳出局部最優解的能力。與以塔板數為優化變量的算法不同,本算法以傳質量作為優化變量,并根據質量衡算反求節點的進出口濃度,可避免迭代計算造成的優化效率低下的問題,顯著提高優化效率。

2.1 初始化

讀取算例的輸入參數,設定網絡的結構參數,包括富/貧流股的分流組數、分流數、節點數。設置種群規模Np,初始狀態下的所有個體均為空結構,并且結構中所有節點對應的變量信息全部歸零。

2.2 連續變量游走

需要進行游走操作的變量包括傳質量、分流比和貧流股流量。依次對所有節點做循環,若節點上存在傳質單元,則該傳質單元的連續變量有γ的概率參與游走。當判定要參與游走時,則分流比游走概率為β,傳質量游走概率為1-β,流量游走概率為δ。分流比、傳質量、貧流股流量分別按照式(15)~式(18)游走。游走公式中1-2r確定隨機的游走正負方向,rΔL確定隨機的游走步長,此公式保證步長是以0為中心的正態分布,即大概率以小步長搜索,小概率以大步長搜索。

經過游走后,若一些傳質單元的傳質量、分流比或流量小于預設的下限Mmin、SPRmin、SPLmin、Lmin,認為這些傳質單元已無存在的必要,因此將這種傳質單元清除,傳質量、分流比、流量以及連接關系全部清零。

2.3 新單元生成

每次迭代有一定概率ε選中任意兩個貧富流股的節點,連接形成一個新的傳質單元,并賦予此傳質單元一個初始的傳質量。若該傳質單元所在的分流上沒有其他傳質單元則需要賦予其初始分流比,若該傳質單元所在的貧流股上沒有其他傳質單元,則還需賦予其初始流量,如式(19)~式(22)所示。

2.4 選擇與變異

如式(23)所示,當生成階段結束后,計算費用F1,并與游走前的費用F0進行比較,選擇費用較低的結構保留至下一次迭代中。然而,為了防止優化陷入局部最優,即便F1<F0,也仍有一定概率ζ接受費用較差的結構,避免算法早熟。

式中,向量Sit表示第it次迭代的網絡結構,包含所有節點的全部變量信息。

3 精細搜索策略的提出

3.1 全局搜索與局部搜索的矛盾

啟發式方法是一類具有貪婪性的隨機搜索方法,這決定了啟發式方法難以兼顧全局搜索和局部搜索。以RWCE 算法優化廢水脫酚算例為例進行說明,設定兩組不同精度的優化參數,對比觀察二者所得結果的特點和差異,流股數據和投資費用計算公式見文獻[12]。保持模型結構參數相同,設定不同的優化參數,參數A 設定較大的步長、進化概率、生成概率以及接受差解概率,而參數B 則將以上參數設定為較小值,經過相同次數的迭代后對比兩套參數所得的結果。當以參數A 進行優化時,在第4.205×108次迭代得到費用為338816 USD·a-1的結果并收斂于此,該結果的結構如圖4 所示。從圖中可知,內部貧流股S1和S2的流量分別為4.9262 kg·s-1和2.2539 kg·s-1,而二者的流量限制分別為5.0000 kg·s-1和3.0000 kg·s-1,由于使用內部貧流股不會帶來費用的增加,因此應盡可能多地使用內部貧流股。由此可知,該結構的連續變量仍有優化的空間,但受限于基礎算法的固定優化參數,難以得到精度更高的結果,算法往往收斂于非局部最優解。因此,出于保持全局搜索能力的目的而設定的固定參數不利于求解更精確的網絡結構,即局部搜索能力不足。

圖4 參數A得到的結構Fig.4 The structure obtained with parameter A

當以參數B 進行優化時,如圖5 所示,得到了一個包含4 個傳質單元的結構。圖6 為兩套參數優化過程的TAC 曲線,可以看出使用參數B 時TAC 下降頻次較少,且下降幅度較小,很快便收斂于368794 USD·a-1,此費用較前者更高。實驗結果說明,全程的高精度的優化對全局優化來說是不利的,不僅優化效率降低,而且一旦陷入局部最優解便難以跳出,導致算法的持續優化能力不足。

圖5 參數B得到的結構Fig.5 The structure obtained with parameter B

圖6 不同優化參數的TAC迭代曲線Fig.6 TAC iteration curves with different optimization parameters

分析可知,產生以上兩個實驗現象的原因如下:首先,MINLP 問題的局部最優解眾多,全局搜索和局部搜索難以兼顧,無論傾向于哪種搜索都會錯過一些更好的解;其次,整個優化過程對精度的要求是不同的,優化初期可以搜索到的可行解數量巨大,對算法全局搜索能力的需求更大,而優化后期很難在全局范圍內搜索到更好的解,但在局部范圍內仍可以搜索到更好的解,因此對算法局部搜索能力的需求更大。鑒于此,可采用前期以較大的步長、進化概率、生成概率,以及接受差解概率進行全局搜索的方式,獲得一個初始解,再通過精度更高的方法對該解做進一步的優化,以彌補初始解由于精度不足而未被充分優化的問題。

3.2 精細搜索策略

提出兩種應用于質量交換網絡優化的精細搜索方法。方法1 的核心思想是:將基礎RWCE 得到的解復制生成大量的新個體,每個新個體均具有相同的初始結構。以該初始結構為中心,在鄰域內以更精細化的隨機游走進行高精度搜索,并將所得到的更優解作為初始解,再次復制給全部個體進行鄰域內的高精度搜索,以此實現更強的局部搜索能力。方法2 的核心思想是:借助確定性方法的高精度和快速收斂的特性,對連續變量進行優化。使用坐標輪換法[22]對基礎RWCE 算法得到的結構進行搜索,該方法將多維的MEN 問題轉化為一維優化問題,并通過黃金分割法求解極小值點,可降低求解難度,同時也擁有較快的求解速度。

精細搜索是以基礎算法獲得的結構為起點的局部搜索,而局部搜索能獲得更好解的前提是基礎算法得到的結構處于較優的解空間內。為了保證基礎RWCE 算法能獲得具有優化潛力的結構,本文通過算法參數的設置來強化初次優化時的全局搜索能力,包括較大的游走步長、傳質單元進化概率、生成概率和接受差解概率,保證算法在全局范圍內獲得一個具有優化潛力的結構,而非收斂于局部最優解,在此基礎之上,再進行精細搜索便可得到更優解。兩種精細搜索方法步驟如下。

方法1

(1) 初始化階段:與初始結構設置相同的模型參數,將結構中所有節點的傳質量、濃度、分流比和流量全部置零。

(2)個體分化階段:讀取基礎RWCE 算法所得最優解的結構數據,并賦值給種群中的所有個體,將該結構的TAC記錄為F0。

(3) 迭代階段:種群中每個個體使用高精度的RWCE 算法依次優化,當迭代次數達到ITf,max時,對比種群最優費用Fp與初始費用F0,如果Fp小于F0,則將種群最優個體的結構Sp記為歷史最優結構SB,否則將初始結構S0記為歷史最優結構,最后將迭代次數歸零。

(4)回代階段:將歷史最優解回代到步驟(2)中,結構數據賦值給種群中的全部個體,作為下一輪迭代的初始結構。

值得注意的是,此方法是基于基礎RWCE 實現的,保留了RWCE 算法的游走、生成、選擇與變異等機制,在以局部搜索為主的同時也并未喪失全局搜索能力,因此不僅可以對連續變量做進一步優化,也具有一定的結構變異能力。此外,方法1 是基于隨機搜索的,收斂速度比較慢,雖然精度較基礎算法有所提高,但得到的解仍不是數學意義上的局部最優,因此在對求解速度和精度有較高要求時,該方法存在一定的不足。操作流程如圖7所示。

圖7 基于啟發式算法的精細搜索策略流程Fig.7 Flow chart of fine search strategy based on heuristic algorithm

方法2

MEN 問題包含傳質量、分流比、流量三種連續變量,它們共同組成n維向量x。本方法對優化問題進行降維處理,依次沿著x的各個分量方向進行一維搜索,同時保持剩余n-1 個分量不變。式(24)~式(26)表示以xk-1為起點,沿著ek分量方向進行一維搜索,求取該方向上的極小值點和最佳步長λB,并將該極小值點作為下一分量的優化起點,對剩下的分量采取相同的操作。當所有分量都優化過一輪之后,再次從第一個分量開始優化,直至某一輪優化前后的目標函數值之差小于η時終止迭代。

若f(α1)<f(α2),則舍棄區間(α2,b],并將區間[a,α2]作為新搜索區間[a1,b1]。若f(α1)>f(α2),則舍棄區間[a,α1),并將區間[α1,b]作為新搜索區間[a1,b1]。重復上述步驟,當搜索區間長度縮短至預先設定長度θ時,取區間中點α*作為極小值點,終止迭代,α*與起始點之差即為最佳步長。

基于啟發式算法的方法1是一種通用的精細搜索方法,可有效處理不同工況的算例,同時保持了一定的結構變異能力。而處理以填料塔為傳質單元的算例時,由于傳質單元投資費用是關于傳質量的連續函數,因此基于確定性方法的方法2 具有較高的精度和較快的收斂速度。但方法2只能對連續變量進行優化,并且只適用于在一個較優解的基礎上做小范圍的調整,不適用于全局優化。

本文求解環境為Intel(R) Xeon(R) Gold 6226R CPU,主頻2.9 GHz,64 GB RAM,采用Fortran90編程。

4 算例驗證

4.1 算例1 焦爐氣脫硫

本算例取自于文獻[3],流股數據如表1所示,主要目的是去除焦爐氣中的硫化氫,以減少燃燒產生的SO2。過程貧流股S1 采用稀氨水作為吸收劑,外部貧流股S2 采用冷凍甲醇作為吸收劑。傳質單元采用篩板塔,投資費用如式(29)~式(34)所示。當相平衡方程為線性時,塔板數的近似計算通常使用Kremser 方程求得,但由于方程形式和對數項的存在,計算過程中可能會產生分母為0的情況,這導致了公式的非連續性。因此本文采用Fraser 等[23]提出的塔板數公式(30),通過Underwood 近似將Kremser方程轉化為連續函數,以此降低計算難度,其中n取1/3,每塊塔板的單價為4552 USD·a-1。貧流股總操作費用由單位操作費用乘以流量求得。RWCE算法參數為γ=0.3,β=0.2,δ=0.01,ε=0.002,ζ=0.01,ΔLSP=0.01,ΔLM=3×10-5,ΔLL=0.01,Mmax=0.002,Lmax=0.002,Mmin=10-5,SPRmin=0.01,SPLmin=0.01,Lmin=0.02。

表1 算例1的流股數據Table 1 Stream data of case 1

使用基礎RWCE 優化此算例,花費時長13278 s得到TAC=408996 USD·a-1的結果,包含4 個傳質單元,總計22 塊塔板,并且S1 流股分為了兩股。此結果已低于現有文獻中的結果,但仍有繼續優化的空間,因此使用精細搜索策略對該結果做進一步優化。采用方法1 經過25441 s 獲得了TAC=407308 USD·a-1的結果,而采用方法2 經過24 s 獲得TAC=407801 USD·a-1的結果,方法1 獲得了費用更低的結果。方法1 的整個優化過程的TAC 曲線如圖8 所示。基礎RWCE 算法前期使TAC 不斷下降,隨后TAC 長期停滯,因此在5.000×108次迭代后使用方法1 進行精細搜索,引入精細搜索策略后TAC 快速下降,并在接下來的優化過程中多次下降。方法1所得結果相比基礎RWCE 算法所得結果TAC 下降了1688 USD·a-1,相比文獻中的最優結果TAC 下降了3257 USD·a-1。精細搜索前后的MEN 結構如圖9所示,相比基礎RWCE 所得結果,S1 的流量下降0.0027 kg·s-1,S2 的流量增加了0.0181 kg·s-1,貧流股操作費用上升了2862 USD·a-1,但傳質單元R1-S1 的塔板數減少了1 塊,節省了4552 USD·a-1的投資費用,因此使TAC 進一步下降。表2 將本文結果與文獻中的結果進行了對比,得到了目前最優的結果,由此可見,精細搜索策略可以實現對已有結構的深度優化,提升了算法的局部搜索能力。

表2 算例1的結果對比Table 2 Result comparison of case 1

圖8 算例1優化全過程的TAC迭代曲線Fig.8 TAC iterative curve of the whole optimization process of case 1

圖9 算例1精細搜索前后的網絡結構Fig.9 Network structure of case 1 before and after fine search

4.2 算例2 空氣除氨

本算例取自于文獻[29],富流股由5 股空氣組成,主要目的是從空氣中去除氨氣。S1、S2 為過程貧流股,S3 為外部貧流股,傳質單元采用填料塔,費用計算方法如式(35)~式(37)所示,流股數據如表3所示。

表3 算例2的流股數據Table 3 Stream data of case 2

式中,618為單位高度投資費用;MASSir,jr,kr表示設備質量;0.225 為年度化系數,年操作時長為8150 h;0.02 為總傳質系數;LMCDir,jr,kr表示對數平均濃度差;Δyleft和Δyright分別表示傳質單元左右兩端的傳質驅動力。RWCE 算法參數為γ=0.3,β=0.2,δ=0.1,ε=0.005,ζ=0.001,ΔLSP=0.03,ΔLM=3×10-4,ΔLL=0.1,Mmax=0.002,Lmax=1,Mmin=5×10-5,SPRmin=0.01,SPLmin=0.01,Lmin=0.02。

高志輝[28]采用自適應模擬退火遺傳算法,得到一個包含2 次分流、9 個傳質單元的結構,年綜合費用為120717 USD·a-1,但從給出的單元流股信息表中發現,5 號單元的富端濃度不符合傳質可行性約束,因此結果存在誤差。本文采用基礎RWCE 算法在1.845×108次迭代時得到TAC=129067 USD·a-1的解,用時21977 s,并且直到5.000×108次迭代也未能實現費用下降,因此使用精細搜索策略對該結構進一步優化。采用方法1進行精細搜索時,經過1.751×108次迭代,用時19263 s,獲得了128065 USD·a-1的結果。采用方法2 進行精細搜索用時51 s,經過1255 次迭代,TAC 下降至127807 USD·a-1,如圖10所示。在本算例中,方法2的速度明顯快于方法1,并且費用也有小幅下降,表現出更好的適用性。本文最優結果相比基礎程序得到的結果下降了1260 USD·a-1,相 比 文 獻[25]中 的 結 果 下 降 了2093 USD·a-1。精細搜索前后的網絡結構如圖11 所示,由于本算例中使用填料塔作為傳質單元,因此圖中圓圈中的數字表示傳質單元的序號而非塔板數。

圖10 算例2基礎RWCE與精細搜索優化過程的TAC曲線Fig.10 TAC curves of basic RWCE and fine search optimization process of case 2

圖11兩個結構間的差異主要在貧流股S2和S3的流量,以及與之連接的傳質單元的傳質量。經過精細搜索后,S2 的流量增加了0.0449 kg·s-1,與之相連的8號單元傳質量增加了0.00058 kg·s-1,與8號單元耦合的9號單元傳質量減少了0.00058 kg·s-1,與9號單元在同一組分流的1號單元傳質量增加了0.00028 kg·s-1,與1號單元耦合的2號單元減少了0.00028 kg·s-1,這導致了外部貧流股S3承擔的傳質負荷降低,因此所需流量也降低,使得操作費用降低了469 USD·a-1。此外,精細搜索后傳質量有了更好的分配方式,使得結構中設備總投資降低790 USD·a-1。與文獻結果的對比列于表4中,驗證了基于確定性方法的精細搜索策略對提升算法局部搜索能力的有效性。

表4 算例2的結果對比Table 4 Result comparison of case 2

圖11 算例2精細搜索前后的結構Fig.11 Network structure of case 2 before and after fine search

5 結 論

(1)低精度的優化參數雖然全局搜索能力較強,但會導致連續變量難以得到充分優化,無法收斂至局部最優解。而高精度的優化參數全局搜索能力較弱,陷入局部最優解陷阱后便難以跳出,導致算法收斂于較差的結果。

(2)提出兩種用于對初始解深度優化的精細搜索策略。方法1將最優結構進行不斷回代和分化,并以高精度的RWCE 算法進行局部搜索,具有適用范圍廣、優化效果好的優點,且保留了一定的全局搜索能力。方法2采用的坐標輪換法對每個變量的最優解逐個求解,而最優解通過黃金分割法求得,具有速度快、精度高的優點,適合使用填料塔的實例。

(3)將精細搜索策略應用于焦爐氣脫硫和空氣除氨算例中,分別得到了407308 USD·a-1和127807 USD·a-1的優化結果,與現有文獻結果相比經濟性提升,驗證了該策略的有效性,為費用最小質量交換網絡問題提供了一種新的行之有效的優化方法。

符 號 說 明

猜你喜歡
優化結構質量
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
“質量”知識鞏固
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
質量守恒定律考什么
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
做夢導致睡眠質量差嗎
論《日出》的結構
主站蜘蛛池模板: 国产最新无码专区在线| 超碰aⅴ人人做人人爽欧美| 五月婷婷精品| 成人午夜亚洲影视在线观看| 亚洲天堂日韩av电影| 丰满人妻久久中文字幕| 三级毛片在线播放| 日韩高清中文字幕| 亚洲第一精品福利| a色毛片免费视频| 久久精品中文字幕少妇| 日本午夜三级| 国产日产欧美精品| 久久五月天综合| 日韩美一区二区| 国产美女精品人人做人人爽| 毛片网站观看| 亚洲第一视频网| 亚洲天堂首页| 玩两个丰满老熟女久久网| 亚洲swag精品自拍一区| 一级福利视频| 97视频精品全国在线观看| 国产精品尹人在线观看| 最新亚洲av女人的天堂| 欧美不卡视频在线观看| 亚洲成A人V欧美综合| 日本黄色不卡视频| 亚洲成人黄色网址| 日韩欧美综合在线制服| 福利在线免费视频| 国产无码精品在线播放| 无码免费视频| 久操中文在线| 九色综合视频网| 国产一区二区三区视频| 99精品国产自在现线观看| 亚洲午夜18| 精品无码一区二区三区电影| 精品国产Av电影无码久久久| 人妻丰满熟妇AV无码区| 久久综合色天堂av| 日本免费一区视频| 热思思久久免费视频| 一本大道无码日韩精品影视| 激情成人综合网| 亚洲视频色图| 亚洲嫩模喷白浆| 九色视频一区| 日韩在线永久免费播放| 伊人激情久久综合中文字幕| 综合社区亚洲熟妇p| 国产超碰在线观看| 亚洲首页国产精品丝袜| 国产伦精品一区二区三区视频优播 | 日韩视频福利| 亚洲国产精品VA在线看黑人| 久久综合国产乱子免费| 亚洲,国产,日韩,综合一区| 女人18一级毛片免费观看| 色婷婷国产精品视频| 中文字幕av一区二区三区欲色| 秋霞午夜国产精品成人片| 国产一二视频| 色婷婷色丁香| 国产拍揄自揄精品视频网站| 999精品色在线观看| YW尤物AV无码国产在线观看| 婷婷六月综合| 国产一在线| 中文字幕在线不卡视频| 亚洲欧美不卡视频| 婷五月综合| 国内精品视频在线| 亚洲成人动漫在线| yjizz视频最新网站在线| 日韩欧美网址| 漂亮人妻被中出中文字幕久久| 成人精品免费视频| 久久99国产精品成人欧美| 日本高清有码人妻| 亚洲综合天堂网|