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

RWCE算法中采用單元重構策略激勵換熱網絡結構優化

2021-06-30 01:29:48韓正恒崔國民趙倩倩肖媛張冠華
化工學報 2021年6期
關鍵詞:優化結構策略

韓正恒,崔國民,趙倩倩,肖媛,張冠華

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

引 言

換熱網絡主要由冷熱物流、換熱器、加熱器、冷卻器等構成,作為進行能量回收利用的重要環節,廣泛存在于化工生產領域。換熱網絡優化屬于混合整數非線性規劃問題,常規求解方法有夾點法[1-2]、確定性方法[3]、啟發式方法[4-6]等,其中基于隨機技術的啟發式方法不僅可以實現整型變量和連續變量同步優化,而且借助計算機編程實現求解過程,操作簡單、易于實現、計算高效,總體表現優于其他方法,受到越來越多學者的青睞。常見啟發式方法有遺傳算法[7]、模擬退火算法[8]、粒子群算法[9]、微分進化算法[10]、強制進化隨機游走算法(random walk algorithm with compulsive evolution,RWCE)[11]等,許多學者將其應用于換熱網絡優化相關研究并取得了豐碩成果[12-14]。

換熱網絡優化涉及連續型變量優化與整型變量優化,常用的啟發式算法可以做到同步優化連續型變量與整型變量。為了獲得更好的優化效果,通常從兩個方面采取措施。一方面是調整算法內部參數或算法執行方式,側重算法內部自身改進:Silva等[15]提出了一種基于粒子群算法的多周期換熱網絡綜合方法;趙亮等[16]建立基于無分流分級超結構模型的雙層優化算法,外層用遺傳算法搜索最優結構,內層用粒子群算法求解連續非線性規劃子問題;陳帥等[17]提出自適應調節速度權重策略改進了粒子群算法;Pav?o等[18]將并行處理技術應用于遺傳算法與粒子群算法,有效提升了其優化效率與優化質量;Xiao等[19]通過調整步長分布尺度拓寬RWCE算法搜索域,增強了算法全局搜索能力;Aguitoni等[20]采用模擬退火算法優化拓撲結構,微分進化算法優化連續變量,優勢互補提升優化質量。另一方面是調整模型或結構攝動,側重算法外部結構改進:Xu等[21]在有分流節點非結構模型中引入串聯結構,豐富結構匹配多樣性;Kim等[22]提出包含流股多重匹配的擴展分級超結構模型,增強分流控制與混合溫度控制,取得其他模型無法得到的結構;韓正恒等[23]提出結構融合競爭策略,通過換熱單元競爭優化的方式挖掘不同個體結構進化潛力;Nair等[24]減少分級超結構模型中的結構約束,允許重復匹配、交叉流、旁路及多個公用工程的存在,可優化出更復雜的換熱網絡結構;Zamora等[25]基于擴展超結構模型,對公用工程重新定位,允許公用工程自由選擇位置,得到了匹配更佳的新結構。上述研究內容中無論是算法內部改進還是外部結構改進均體現為換熱網絡結構進化,相較而言對結構優化過程直接處理的外部結構改進是更直觀的處理方法。

RWCE算法優化換熱網絡時展現出較強的性能,但仍具有一定的改進空間。現有的研究基于上述算法和結構兩個出發點,主要側重根據算法后期表現提出改進方法,以更合理的換熱網絡結構為優化目標,卻缺乏對優化過程中算法在結構優化上的作用效果進行較為具體的研究。本文從該角度出發,通過監測換熱單元的進化過程,分析算法優化的作用特點,指出影響結構優化的主要因素。提出換熱單元重構策略改進算法優化流程,以較直接的方式突破結構桎梏,并建立了進化狀態實時監測指標,指導策略執行,求解費用更低的換熱網絡結構。

1 換熱網絡優化模型

采用節點非結構(node-wise non-structural model,NW-NSM)模型[26]優化換熱網絡,該模型初期為空結構,僅在流股上預設一定數量的節點用于匹配換熱單元。如圖1所示,該換熱網絡包含2條熱流股和2條冷流股,每條流股上6個節點,優化過程中在熱流股和冷流股上分別隨機選擇一個節點匹配形成換熱單元,若該單元滿足約束條件則認為匹配成功。圖1中包含4個換熱單元,流股末端分別為冷公用工程與熱公用工程。

圖1 節點非結構模型示意圖Fig.1 NW-NSMdiagram

以最小年綜合費用(total annual cost,TAC)為優化目標,在式(1)中用F表示,包含換熱單元固定投資費用、面積費用及公用工程費用。

式中,LH和LC分別表示熱、冷流股數目。

構造換熱網絡時遵循傳熱學定律及原理,需滿足相應的約束,如傳熱平衡約束、最小溫差約束等。基于各種約束條件進行的目標函數計算決定了換熱網絡優化是混合整數非線性規劃問題,具體條件及計算方法見文獻[26]。

2 換熱網絡結構優化進程分析

2.1 強制進化隨機游走算法

RWCE算法作為一種新型啟發式方法,個體獨立進化,不依賴種群信息交流,且具有獨特的接受差解機制,與其他常用算法相比,具有魯棒性更好、全局搜索能力更強的特點。該算法與節點非結構模型契合度較高,將其應用于節點非結構模型可以充分發揮其搜索能力。基于節點非結構模型的RWCE算法在優化換熱網絡的過程中,通過隨機選擇流股節點匹配形成換熱單元,并基于既有最優結構進行再優化,同時配合接受差解機制,實現廣泛搜索域內的換熱網絡優化。優化流程主要包括換熱單元的生成與進化、個體選擇與變異等步驟,其流程圖如圖2所示。由于節點非結構模型初期為空結構,換熱單元的生成位置完全隨機,理論上可匹配形成任何結構,需要對換熱單元生成與進化過程進行跟蹤探究,以便在更深層次上掌握該模型中換熱網絡結構優化的實際特點。

圖2 基于節點非結構模型的RWCE算法流程圖Fig.2 Flow chart of RWCEbased on NW-NSM

2.2 換熱單元生成與進化特點

換熱網絡結構優化以換熱單元的生成與進化作為基礎,為了探究結構優化特性,應用RWCE算法優化多個算例,通過考察換熱單元在優化進程中的熱負荷變化反映結構優化時的優勢與不足。以15SP1算例[27]為例,該算例中共有8條熱流股、7條冷流股,選取最優個體結構中4個具有代表性的換熱單元,優化過程中其熱負荷變化如圖3所示,500萬步后優化陷入停滯,直至1200萬步最優個體結構與換熱單元熱負荷均不再變化,因此圖中僅顯示前600萬步的變化情況。

圖3 最優結構換熱單元熱負荷變化Fig.3 Heat load variation of heat exchange units in optimal structure

由圖3可以看出,前期優化過程中結構波動較大,存在換熱單元快速生成與消去的過程:1號與2號換熱單元在生成后經過進化迅速穩固地位,在結構中長期存在;3號與4號換熱單元在較短時期內經歷了生成、進化與消去的過程。在100萬~400萬步的優化區間內,1號與2號換熱單元負荷出現小幅度優化。約450萬步時結構再次出現明顯改變,其中,2號換熱單元迅速消去,3號換熱單元再次生成并迅速進化穩定存在,1號換熱單元再度進化,熱負荷增加,形成更穩定的換熱單元。綜合觀察整個優化進程,在前期部分換熱單元迅速占據優勢后,易造成先入為主的現象。上述最優結構中共12個換熱單元,其中8個換熱單元是在進化早期形成的換熱單元,貫穿整個優化進程,經歷持續進化卻一直穩定存在。而期望出現的結構更新與變異過程在約450萬步時出現一次,說明優化后期RWCE算法具有促使結構進一步優化的能力,但概率較小造成中間經歷了較長時間的結構進化停滯期。所以總體優化特點是前期優化效率較高,但部分換熱單元先入為主在結構中占據重要地位,造成中后期結構更新與進化能力不足,最后進入優化停滯階段。因此依靠算法本身的優化方法難以維持較強的結構進化活性,需要借助一定的改進手段實現結構變異,突破先入為主換熱單元的“封鎖”,提升優化質量。

3 換熱單元重構策略

基于個體當前的結構進化特點,對先入為主的換熱單元進行處理是改變結構布局的可行方法。下文提出一種攝動方法消去結構中部分既有換熱單元,為其他換熱單元提供進化空間,對作用效果進行初步探究并建立結構進化狀態評價指標,以合理調控優化時機,發揮改進策略的優勢。

3.1 換熱單元攝動

采用RWCE算法優化15SP1算例,并對換熱單元進行定期攝動,具體操作是每隔50萬步,消去結構中3個換熱單元,目的在于對現有結構形成一定程度的破壞,再利用算法自行生成若干換熱單元,構成新結構并考察其實際進化情況。以年綜合費用下降情況反映加入攝動操作前后的作用效果,如圖4所示。根據圖4,加入攝動操作后,50萬步時,結構攝動產生影響,與正常優化流程相比,一定程度上拖慢了優化進程;100萬步時,結構中的攝動作用帶來明顯的增益效果,約120萬步時,費用明顯優于原優化流程中的同時期費用;后續在200萬與350萬步時,結構攝動作用再次促進了結構優化。同樣優化至600萬步時,正常優化流程所得最優結構如圖5所示,年綜合費用為1593975 USD·a-1;加入攝動操作后所得最優結構如圖6所示,年綜合費用為1537891 USD·a-1。圖中藍色換熱單元為兩結構中具有差異的換熱單元,經過對比可以看出,攝動操作通過改善結構布局實現結構變異與進化,是提升優化效果的有效措施。

圖4 加入攝動前后年綜合費用變化曲線Fig.4 TACcurves before and after the perturbation

圖5 加入攝動前結構圖(1593975 USD·a-1)Fig.5 Structure diagrambefore perturbation

圖6 加入攝動后結構圖(1537891 USD·a-1)Fig.6 Structure diagram after perturbation

3.2 換熱單元重構策略

RWCE算法優化換熱網絡的過程中個體獨立進化,互相之間不進行信息交流。常規優化過程中一般會設置一定的種群規模,主要目的是搜索得到不同個體中的多樣化結構,尋求更全面的優化。為了得到更好的結果,往往需要設置較大的種群規模,優化效率較低。由2.2節的研究可知,先入為主的換熱單元決定了個體主要結構,長期優化過程中,種群規模大并不能對個體進化過程產生影響,反而造成大量無效計算。3.1節進行的結構攝動操作對陷入優化停滯的結構進行有效改進,增強了個體結構變異能力,求解出了年綜合費用更低的新結構。這個過程中實質上是通過提升個體的結構多樣性間接提升了種群多樣性,即單獨個體結構經過攝動可能變異為其他個體的結構,取代了部分種群規模的作用。

綜上所述,提出換熱單元重構策略改進RWCE算法(random walk algorithm with compulsive evolution with unit-reconfiguration strategy,UR-RWCE),借助約束保障較強的結構變異能力,策略執行步驟如下。

步驟一:消去換熱單元。一定周期、一定概率τ0對個體結構進行處理,選定結構中某換熱單元,若隨機數τ<τ0,則消去該換熱單元,熱負荷清零并解除節點匹配關系。若結構中共有NS個換熱單元,則一次處理過程中約消去(NSτ0)個換熱單元。處理方式:

式中,Qdh為被選中消去的換熱單元熱負荷,kW;MCdh為熱流股第dh個節點所連接的冷流股節點編號;MHdc為冷流股第dc個節點所連接的熱流股節點編號;τ0為換熱單元重構概率;τ為0~1之間均勻分布的隨機數。

步驟二:重構換熱單元。若依據式(2)消去了某換熱單元,則隨機選擇一個熱流股節點和一個冷流股節點構建一個新換熱單元,若該換熱單元的匹配關系和溫位均與被消去的換熱單元相同,則放棄在該處構建換熱單元,重新選擇節點構建換熱單元,新換熱單元熱負荷依據原換熱單元熱負荷進行隨機賦值。執行方式如式(3)所示。

式中,Qdh'為新換熱單元熱負荷,kW;Qdh為被消去的換熱單元熱負荷,kW;NH'、NC'為新換熱單元所連接的熱、冷流股;NH、NC為原換熱單元所連接的熱、冷流股編號;T′W—為新換熱單元的溫位;TW為原換熱單元的溫位;λ為0~1之間均勻分布的隨機數。

3.3 換熱網絡結構進化狀態評價指標

結構變異后通常需要一定的時間進行充分優化,若僅采用定期策略攝動的辦法,易出現結構仍處在良好進化狀態卻被再一次攝動改變的情況。對結構進化狀態進行實時監測可以更準確地把握策略攝動時機,當監測到結構長時間未產生有效變異進化時,再對其進行攝動處理。基于三參數Logistics模型建立個體進化狀態實時監測指標:

式中,S(t)為結構進化停滯程度;t為進行結構優化的累積計算步;a決定曲線的區分度,影響S(t)到達臨界值的速度;b決定計算步的累積周期;c決定S(t)的下限。取c=0,a=0.02,b=300時,曲線走勢如圖7所示。當S(t)值略大于0時,代表結構產生有效進化,需要一定時間進行充分優化;當S(t)值趨近1時,表示結構已經進行了充分優化,進入穩定階段,亟須進行結構變異優化。

圖7 三參數Logistics模型變化曲線Fig.7 Three parameter Logistics model curve

設置小于1的臨界值ζ,當S(t)值到達臨界值時,則執行策略,此時t值即為策略執行周期。a、b、ζ賦予不同值時,S(t)到達臨界值ζ的步數是不同的,基于此可以實現對整體進化過程調控,同一個體不同階段的擾動周期、不同個體的擾動周期均不同,增強了結構變異進程的多樣性。t初值為0,根據式(3)進行取值累積計算,即當費用下降或策略執行后,t值歸0,否則t值逐步累加。根據上文研究,當結構費用經過50萬步沒有下降時,一般認為該結構已進入優化停滯狀態,a、b、ζ根據式(6)~式(8)取值,周期主要控制在40萬~60萬步,保障個體優化進程具有不同頻率的攝動。改進策略執行流程如圖8所示。

圖8 策略執行流程圖Fig.8 Flow chart of strategy implementation

4 算例驗證

4.1 算例1

算例1采用20SP算例,13股熱流、7股冷流,該算例取自文獻[18]。該算例流股熱容流率差異較大,流股匹配困難,Pav?o等多位學者應用該算例進行了換熱網絡優化研究[18,28-31]。Zhang等[29]應用改進的布谷鳥搜索算法(CS)取得TAC為1418981 USD·a-1的結構;Xu等[30]采用固定投資松弛策略改進RWCE算法(RSFCC-RWCE),優化該算例取得TAC為1412801 USD·a-1的結構;Rathjens等[31]應用高效的遺傳算法(GA)優化組合策略取得TAC為1407203 USD·a-1的結構;陳子禾等[32]基于多個體平行搜索思想進行參數輪換尋優(DC-RWCE),得到TAC為1401958 USD·a-1的換熱網絡結構,是目前已發表文獻中的最優結果。應用策略改進的RWCE算法(UR-RWCE)優化該算例,設置進化步長L=100 kW,新換熱單元熱負荷初值Q0=200 kW,最小熱負荷限制值Qmin=10 kW,接受差解概率δ=0.01。換熱單元重構概率τ0分別取0、0.2、0.4、0.6、0.8、1,優化結果對比如表1所示。其中重構概率0即RWCE基礎算法優化所得結果,變異比指新結構與基礎算法所得結構對比,具有結構差異的換熱單元數與總換熱單元數的比值。

表1 不同重構概率的優化結果對比Table 1 Comparison of results with different reconfiguration probabilities

根據表1數據,重構概率不宜過大或過小,取較為適宜的中間值,既對結構起到一定擾動作用,又不會造成結構擾動過度。表1中重構概率為0.4時所得結構年綜合費用最低,而年綜合費用與變異比無直接關聯關系。在表1所得結果的基礎上將重構概率設為0.3和0.5,所得結構年綜合費用分別為1407226、1406654 USD·a-1,所以重構概率為0.4所得的結構更優。基礎算法所得結構如圖9所示,重構概率0.4所得結構如圖10所示。圖10中藍色換熱單元表示有效變異的換熱單元,對比可以發現,結構變異部分既包含匹配關系的調整,也包含溫位的調整,促進了流股的合理匹配與熱負荷的合理分配。新結構年綜合費用與原最優結構費用相比有了大幅下降,驗證了改進策略的有效性。本文優化結果與部分文獻較優結果的對比如表2所示。

表2 算例1優化結果與文獻值的比較Table 2 Comparison of optimal resultsfor case 1

圖9 1428046 USD·a-1對應的換熱網絡結構Fig.9 HENstructure with cost of 1428046 USD·a-1

圖10 1395724 USD·a-1對應的換熱網絡結構Fig.10 HENstructure with cost of 1395724 USD·a-1

4.2 算例2

算例2采用15SP2算例,含10股熱流、5股冷流,取自文獻[33]。曹美等[33]探究換熱網絡結構優化過程中的交叉結構并采用禁忌策略改進算法(CSPRWCE),優化該算例取得TAC為5233287 USD·a-1的結構;趙倩倩等[34]分析了優化過程中新換熱單元的熱負荷與生成概率對優化效果的影響,提出參數協同取值策略(CSDP-RWCE)實現了算法優化能力的提升,優化該算例取得TAC為5177785 USD·a-1的結構;趙倩倩等[35]通過調整最大游走步長與換熱單元最小熱負荷約束,改進RWCE算法優化進程(PRRWCE),取得TAC為5169883 USD·a-1的換熱網絡結構,是當前已發表文獻的最好結果,其結構如圖11所示。為了進一步證明本文改進策略的有效性,針對文獻結構繼續優化,重構概率為0.3時得到TAC為5113717 USD·a-1的新結構,如圖12。

圖11 5169883 USD·a-1對應的換熱網絡結構Fig.11 HENstructure with cost of 5169883 USD·a-1

在策略執行優化過程中進行跟蹤監測,結合優化結果進行分析,發現選擇合適的換熱單元進行重構非常重要。對比較好結果和較差結果的優化過程,一般選擇重構換熱溫位匹配不合理的換熱單元更有利于結構進化,因為溫位匹配合適(即高溫位匹配高溫位,低溫位匹配低溫位)有利于降低換熱單元面積。而重構概率的取值一般不宜過大或過小,過小時對結構攝動作用不足,過大時對結構破壞過度,結構中原有較優搭配也被破壞,引起結構變差,因此,再次驗證重構概率適宜取偏中間值。將算例優化前后的結構對比,發現部分換熱單元在重構后順利保留下來,如圖12中藍色換熱單元;部分換熱單元經過重構所得換熱單元在經歷一定優化過程后被消去,依靠算法自身的優化能力生成較為穩定的新換熱單元,如圖12中紅色換熱單元;部分換熱單元未經歷重構過程,如圖中黑色換熱單元,占比較小。藍色換熱單元與原結構匹配關系與溫位不同;紅色換熱單元的出現,說明部分重構換熱單元并不合理,造成結構費用較高,進而被新生的更合理換熱單元取代,但換熱單元的重構過程促進了這部分換熱單元的形成。少部分紅色換熱單元與原結構中部分換熱單元匹配關系和溫位均相同,說明存在一些流股匹配關系在多種結構中均屬于較合理的匹配關系。綜合比較之下,換熱單元重構策略可以有效提升結構多樣性。該算例優化結果與文獻結果對比如表3所示。

圖12 5113717 USD·a-1對應的換熱網絡結構Fig.12 HEN structure with cost of 5113717 USD·a-1

表3 算例2優化結果與文獻值的比較Table 3 Comparison of optimal results for case 2

4.3 算例3

基于上述研究與分析,進一步采用其他算例驗證改進策略的有效性。算例3為9SP算例,4條熱流股、5條冷流股,取自文獻[36]。該算例應用同樣較為廣泛,Xu等[30]提出的固定投資松弛策略優化該算例得到TAC為2917682 USD·a-1的結構,是目前文獻中的最優結果。基礎RWCE算法優化該算例得到TAC為2928507 USD·a-1的換熱網絡結構,采用換熱單元重構策略并設置重構概率為0.3時,得到TAC為2909121 USD·a-1的結構,相比文獻最優結果下降了8561 USD·a-1,再次驗證了換熱單元重構策略的有效性。最優結構如圖13所示,優化結果與文獻結果的對比如表4所示。

表4 算例3優化結果與文獻值的比較Table 4 Comparison of optimal resultsfor case 3

圖13 2909121 USD·a-1對應的換熱網絡結構Fig.13 HENstructure with cost of 2909121 USD·a-1

4.4 算例4

算例4為15SP3算例,含8條熱流股、7條冷流股,該算例與2.2節的15SP1算例流股數量相同,但參數不同,取自文獻[20]。該算例是一個比較新的算例,Aguitoni等[20]應用雙層優化算法(SA-DE)優化該算例所得結構TAC為12389890 USD·a-1。應用RWCE算法優化該算例得到TAC為12060224 USD·a-1的結構。RWCE算法結合應用換熱單元重構策略,設置重構概率0.4時得到TAC為11549145 USD·a-1的換熱網絡結構,如圖14。與文獻結果相比,新結構的TAC值下降了840745 USD·a-1,優化結果與文獻結果的對比如表5所示。通過以上多個算例的驗證,采用換熱網絡重構概率的RWCE算法結構優化能力明顯強于基礎算法,且可以取得目前文獻中的最優結果,充分證明了改進策略的有效性與普適性。

圖14 11549145 USD·a-1對應的換熱網絡結構Fig.14 HENstructure with cost of 11549145 USD·a-1

表5 算例4優化結果與文獻值的比較Table 5 Comparison of optimal resultsfor case 4

5 結 論

(1)RWCE優化換熱網絡結構的特點是前期優化效率較高,結構進化幅度較大,但易出現部分換熱單元在整個優化進程中長期存在并占據穩固地位,形成先入為主的現象,中后期結構更新變異難度較大,長期處在優化停滯狀態。

(2)換熱單元重構策略可以有效促進結構進化,主要作用方式為改變部分換熱單元匹配關系或溫位。但并非所有重構換熱單元都是合理的,不合理的換熱單元會被基礎RWCE算法新生的換熱單元取代。在重構策略的作用下,整體結構經過變異與再進化,形成年綜合費用更低的新結構。

(3)與基礎算法所得結果及文獻結果相比,采用單元重構策略優化實際算例具有明顯優勢。優化算例9SP、20SP、15SP2、15SP3所得結果與文獻最優結果相比均有不同程度的降低,是目前文獻中的最優結果,驗證了本文提出的改進策略應用于換熱網絡結構優化問題的優越性。

符號說明

A——換熱面積,m2

CA——面積費用系數,USD·m-2·a-1

CF——換熱單元固定投資費用,USD·a-1

CHU,CCU——熱、冷公用工程單位能耗費用,USD·kW-1·a-1

MC——與熱流股節點連接的冷流股節點

MH——與冷流股節點連接的熱流股節點

PH、PC——每條熱、冷流股上節點數

Q——換熱單元熱負荷,kW

r1,r2,r3,r4,r5——0—1之間均勻分布隨機數

S——結構進化停滯程度

Z——值為1或0的整型變量

ε——面積費用指數

下角標

i——熱流編號

j——冷流編號

猜你喜歡
優化結構策略
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
例談未知角三角函數值的求解策略
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
我說你做講策略
高中數學復習的具體策略
數學大世界(2018年1期)2018-04-12 05:39:14
論《日出》的結構
主站蜘蛛池模板: 天天躁夜夜躁狠狠躁图片| 天天躁夜夜躁狠狠躁图片| 免费中文字幕在在线不卡| 亚洲国产精品VA在线看黑人| 黄色一级视频欧美| 精品久久蜜桃| 免费国产不卡午夜福在线观看| 91啪在线| 91娇喘视频| 日本在线视频免费| 日韩二区三区无| 欧美劲爆第一页| 国产人成在线观看| 婷婷午夜影院| 91尤物国产尤物福利在线| 国产午夜小视频| 韩国v欧美v亚洲v日本v| a天堂视频| 午夜精品区| 久久99精品久久久大学生| 欧美成人国产| 婷婷色丁香综合激情| 亚洲国产亚洲综合在线尤物| 国产精品一区二区在线播放| 国产主播福利在线观看| 丁香五月激情图片| 99久久精品免费看国产免费软件 | 一级毛片在线播放免费| 特级做a爰片毛片免费69| 黄色在线不卡| 一本视频精品中文字幕| 无码免费试看| 五月天福利视频 | 精品久久久无码专区中文字幕| 无码AV日韩一二三区| 欧美日韩第二页| 手机成人午夜在线视频| 波多野结衣视频一区二区| 久久久久久尹人网香蕉| 91成人在线观看| 欧美中出一区二区| 国产日韩久久久久无码精品| 欧美成人一级| 亚洲一区二区三区麻豆| 一级爆乳无码av| 国产高清在线观看91精品| 高清无码一本到东京热| 国产真实乱子伦视频播放| av一区二区三区在线观看| 国产成人免费观看在线视频| 欧美成a人片在线观看| 欧美一区二区丝袜高跟鞋| 精品视频一区二区观看| 欧美亚洲欧美区| 怡春院欧美一区二区三区免费| 97se亚洲| 99视频精品全国免费品| 91无码视频在线观看| 久久天天躁狠狠躁夜夜躁| 亚洲第一视频网站| 国产人成网线在线播放va| 久久国产拍爱| 国产精品区网红主播在线观看| 美女高潮全身流白浆福利区| 91高清在线视频| 91无码人妻精品一区二区蜜桃| 亚洲美女久久| 亚洲国产中文在线二区三区免| 人妻精品久久无码区| 国产欧美日韩综合在线第一| 网友自拍视频精品区| 日本免费福利视频| 欧美成人影院亚洲综合图| 激情午夜婷婷| 久久综合伊人77777| 曰AV在线无码| 99re热精品视频中文字幕不卡| 婷婷色婷婷| 最新亚洲人成网站在线观看| 精品视频在线一区| 亚洲不卡影院| 欧美性猛交一区二区三区|