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

板翅式換熱器鋸齒型翅片參數的遺傳算法優化研究

2015-03-07 02:14:54楊輝著文鍵童欣李科厲彥忠王斯民
西安交通大學學報 2015年12期
關鍵詞:關聯優化

楊輝著,文鍵,童欣,李科,厲彥忠,王斯民

(1.西安交通大學能源與動力工程學院, 710049, 西安;2.西安交通大學化學工程與技術學院, 710049, 西安)

?

板翅式換熱器鋸齒型翅片參數的遺傳算法優化研究

楊輝著1,文鍵1,童欣1,李科1,厲彥忠1,王斯民2

(1.西安交通大學能源與動力工程學院, 710049, 西安;2.西安交通大學化學工程與技術學院, 710049, 西安)

針對目前鋸齒型板翅式換熱器未能同時優化多參數,或者大多優化研究存在對經驗關聯式依賴的問題,提出了利用Kriging響應面來近似目標函數與設計變量之間的關系、應用遺傳算法對鋸齒型板翅式換熱器翅片結構參數的優化方法。在維持翅片通道雷諾數為800時,把換熱器的最大j因子、最小f因子和最大FTEF因子作為3個單目標函數,對翅片的翅片高度h、翅片間距s、翅片厚度t和翅片節距l進行了優化研究。研究結果表明:翅片高度h與翅片間距s對換熱器綜合性能FTEF因子呈正增長,而翅片厚度t和翅片節距l呈負增長;在翅片高度為9.5 mm、翅片間距為2.2 mm、翅片厚度為0.1 mm和翅片節距為3 mm時,換熱器性能最佳;結合Kriging響應面的遺傳算法克服了傳統優化方法對經驗關聯式的依賴。該研究結果可以指導鋸齒型板翅式換熱器的優化設計。

板翅式換熱器;遺傳算法;Kriging響應面;優化設計

板翅式換熱器的基本結構是,在相鄰兩隔板之間放置的二次傳熱面構成狹窄的矩形流體通道,其傳熱與流動阻力性能主要決定于具有二次傳熱的翅片表面特性。作為緊湊式換熱器的基本傳熱元件,此二次傳熱面可以擴大傳熱面積,提高傳熱效率,增強結構緊湊性,提升換熱器強度及承壓能力。根據不同的工藝,翅片的型式主要有平直翅片、多孔翅片、鋸齒翅片、波紋翅片和百葉窗翅片[1]。其中,鋸齒型翅片以其緊湊、高的可靠性及傳熱效率高等優點得到了廣泛的應用,不過也存在流動阻力大的缺點。因此,對鋸齒型板翅式換熱器進行優化研究,對于提高其傳熱量和降低運行費用具有重大的意義。

目前,板翅式換熱器的設計主要依據經驗選擇、多次驗算和前人所得結果,通過簡單的流動與傳熱理論計算,在滿足設定的換熱性能和阻力要求時,即可作為所參考的換熱器形式和翅片類型,然而這樣得到的換熱器形式并不能保證是最佳的方案。遺傳算法[2]的基本理論依據是模擬生物在自然環境中的遺傳和進化過程而形成的一種自適應全局搜索算法,在板翅式換熱器優化研究中得到廣泛的應用。文獻[3-5]以板翅式換熱器重量和效率為目標函數,文獻[6]以板翅式換熱器火積耗散最小為目標函數,文獻[7]以板翅式換熱器熵產數為目標,對換熱器進行了單目標遺傳算法的優化研究。此外,文獻[8]以能量消耗與材料花費為目標函數,文獻[9]以換熱量與壓力損失為目標函數,文獻[10]以總換熱量與系統全年運行總成本為目標函數,對板翅式換熱器進行了多目標遺傳算法的優化研究。然而,上述的優化方法基本上都需要依靠經典關聯式計算翅片的傳熱與阻力性能。文獻[11-19]通過實驗或數值研究的方法,提出了鋸齒型翅片表面特性的關聯式,但是這些關聯式的工質都是空氣,將其應用于大Pr為工質的情況,存在很大誤差。目前,鋸齒型板翅式換熱器還嚴重缺乏以大Pr為工質的經驗關聯式。

本文以鋸齒型板翅式換熱器為研究對象,首先以翅片高度、翅片間距、翅片厚度和翅片節距為優化參數進行數值研究,然后結合Kriging響應面技術,基于遺傳算法,在維持翅片通道Re恒定時,以j因子、f因子和FTEF因子為目標函數,對換熱器的翅片高度、翅片間距、翅片厚度和翅片節距進行了優化設計。

1 計算模型與數值方法

1.1 優化設計物理模型

圖1為鋸齒型板翅式換熱器翅片結構參數軸測圖,其中翅片參數包括翅片高度h、翅片間距s、翅片厚度t和翅片節距l。圖2為換熱器翅片的平面圖,可知單周期翅片端部參與換熱的面積為Aend=2t(h-t)+t(s-2t),所以本文定義的當量直徑[19]為

(1)

針對目前鋸齒型板翅式換熱器主要使用的翅片參數范圍,此次選取的翅片優化范圍是:h為4.7~9.5 mm,s為1.4~3 mm,t為0.1~0.5 mm,l為3~9 mm。

圖1 鋸齒型板翅式換熱器翅片結構參數的軸測圖

圖2 鋸齒型板翅式換熱器翅片的平面圖

1.2 基本方程及數值方法

本文的基本方程包括質量、動量以及能量方程。連續性方程為

(2)

動量方程為

(3)

能量方程為

(4)

模擬計算時,入口為質量流量邊界條件,其對應的翅片通道Re維持為800,入口溫度為300 K;出口為壓力出口;假定上下隔板被充足的飽和蒸汽加熱,即設定隔板表面為定壁溫(373.15 K);側面為周期邊界條件;流固接觸面為默認的耦合邊界條件。板翅式換熱器材料為鋁,且忽略與外界的輻射與對流傳熱。通道流體為空氣,其物性假定為常物性。壓力-速度耦合采用SIMPLE算法,動量及能量方程使用二階迎風格式。能量方程的計算殘差為1×10-10,其他參數的殘差為1×10-6。

圖3 鋸齒型翅片網格示意圖

計算模型的網格由mesh模塊生產,全部采用結構化網格,同時對壁面劃分了邊界層網格及局部加密處理,板翅式換熱器網格如圖3所示。同時,為了保證計算結果的準確,對生成的網格進行了獨立性驗證。圖4為對應無量綱y+為1、1.5、2、3、4和5的6組71JC2203翅片的網格無關性驗證圖。由圖可知,在網格數增加到3 703 700后,換熱器j因子和f因子變化都小于1%,考慮到計算的精度和縮短計算時間,選擇網格為3 703 700。本文所有的鋸齒型板翅式換熱器的網格單元數為2 410 260~4 102 134。

圖4 71JC2203翅片的網格無關性驗證

1.3 數值模擬有效性的驗證

為了驗證數值模擬的正確性,將模擬計算值與文獻[10]鋸齒型翅片的實驗值進行對比,結果如圖5所示,其中翅片型號為95JC1803。由圖可知,j因子和f因子的數值結果與實驗數據吻合較好。j因子數值結果與實驗值的均方根誤差為8.56%,而f因子數值結果與實驗值的均方根誤差為7.91%。由文獻[10]得知,j因子及f因子存在5%的實驗系統誤差,因此可認為數值模擬結果是正確的。

圖5 數值模擬與實驗結果的對比

2 優化研究

2.1 優化原理

圖6為遺傳算法的優化過程。優化過程包括初始試驗設計、構造響應面和優化目標。在試驗設計的選擇中采用具有試驗次數少、精度高、預測性好的中心組合設計。基于統計設計的先進優化方式,善于體現影響因素之間的交互作用程度,可以評估因素的非線性影響。響應面方法是指構造源模型的響應面,以解決源模型的設計或分析等問題的近似方法。這種利用響應面來近似目標函數與設計變量之間的關系,代替源函數進行估值,可以大大減少源函數的估值次數[20],克服了傳統的全局優化方法對源函數的估值次數多、效率低的問題。Kriging方法是一種通過已知點來預測未知觀察點的插值方法,其利用方差的變化來表達空間的變化,可以保證由空間分布得到的預測值的誤差最小[21]。Kriging差值方法表示為Z(x)=f(x)+g(x),其中f(x)為二階多項式方程(決定了模型了的全局行為),g(x)為擾動項(決定了模型的局部行為)。Kriging響應面可以同時提供目標點處的預測值和預測誤差,保證響應面通過所有的樣本點,使全局優化算法具有良好的全局性和局部性。最后,利用Kriging構造的響應面,通過遺傳算法優化目標函數。這種基于響應面技術的遺傳算法的優化研究,是利用響應面創建的輸入、輸出函數,而非傳統的經驗關聯式,所以在處理無經驗關聯式,或者經驗關聯式不夠準確的情況非常有效。不過這種方法得到的優化解是近似的,需要對優化解進行驗證。

圖6 遺傳算法的優化過程

2.2 數據處理

計算Nu及j因子

(5)

(6)

(7)

(8)

(9)

式中:u為通道速度;D為當量直徑;h為傳熱系數;m為流體質量流量;cp為比定壓熱容;T為溫度;Δtm為對數平均溫差;λ為導熱系數;μ為動力黏度;下標in、out分別表示進、出口。

計算摩擦因子f

(10)

式中:Δp為壓差;ρ為流質密度;L為翅片長度。

計算FTEF因子

(11)

3 結果分析與討論

3.1 數值結果與經驗關聯式的對比

(a)j因子

采用中心組合設計的方法,均勻分配了50個設計點。為了驗證數值結果的正確性,將本文結果與Mangles & Bergles關聯式和Kim & Lee關聯式進行了對比,結果如圖7所示,其中橫坐標為關聯式的計算結果j1、f1,縱坐標為本文結果。由圖7a可知,Mangles & Bergles和Kim & Lee關聯式的j因子與數值結果吻合的較好。Mangles & Bergles關聯式的誤差范圍在±15%,均方根誤差為6.0%,而Kim & Lee關聯式的誤差范圍在±10%,均方根誤差為4.7%。由圖7b可知,Mangles & Bergles和Kim & Lee關聯式的f因子與本文結果的誤差基本維持在±20%之間,不過也存在少量Mangles & Bergles關聯式的誤差達到25%。Mangles & Bergles關聯式與本文結果的均方根誤差為18.7%,而Kim & Lee關聯式的均方根誤差為15.4%。產生誤差的原因可能是f因子與翅片當量直徑成正比,j因子與當量直徑無關。Mangles & Bergles和Kim & Lee關聯式定義的當量直徑比本文的大,因此本文結果的f因子必然比經驗關聯式的值小。此外,數值模擬忽略翅片表面的粗糙度,也導致f因子偏小。由于經驗關聯式本身也存在±20%誤差,可認為本文數值模擬結果是準確的。

(b)f因子圖7 本文結果與經驗關聯式結果的對比

3.2 翅片參數對換熱性能的影響

3.2.1 翅片高度h的影響
圖8為s=2.2 mm、t=0.3 mm、l=6 mm時,不同翅片高度h對j因子、f因子及FTEF因子的影響。由圖8可知,隨著h的增加,j因子、f因子及FTEF因子都逐漸增加,其中j因子增加23.8%,f因子增加14.3%,FTEF因子增加18.4%。隨著h增加,翅片二次傳熱面積顯著增加,增大了流體與翅片的換熱面積,強化了傳熱。雖然h增加使翅片通道流速減小,但是在維持通道Re一定時,流道內的質量流量是增加的,最終導致j因子、f因子都所有增加。增加h有利于強化換熱,所以在滿足允許的壓降下,應適當增加翅片高度。

圖8 翅片高度h對換熱性能的影響

圖9 翅片間距s對換熱性能的影響

3.2.2 翅片間距s的影響
圖9為h=7.1 mm、t=0.3 mm、l=6 mm時,不同s對j因子、f因子及FTEF因子的影響。由圖9可知,隨著s的增加,j因子和FTEF因子先有所增加后趨于穩定,而f因子先快速減小后變緩,其中j因子增加16%,f因子降低47%,FTEF因子增加27.9%。隨著s增加,翅片一次表面和通道內的質量流量都增加,強化了傳熱,使j因子增加。翅片通道內流速快速下降使f因子顯著下降,同時通道內流量增加又減緩了f因子的下降。低翅片間距使通道壓降變大,翅片性能變差,應該適當設計大的翅片間距。

3.2.3 翅片厚度t的影響
圖10為h=7.1 mm、s=2.2 mm、l=6 mm時,不同t對j因子、f因子及FTEF因子的影響。由圖10可知,隨著t的增加,j因子幾乎不變,FTEF因子逐漸減小,而f因子迅速增加,其中,FTEF因子減小27.8%,f因子增加88.6%。翅片厚度增加導致通道流通面積減小,流速快速增加,引起f因子的顯著提高。同時,雖然流速增加有利于換熱,但是通道質量流量和換熱面積都減小,結果j因子幾乎不變。翅片厚度對換熱器的壓降具有特別明顯的影響,在低壓板翅式換熱器的設計中應該選擇薄翅片,而高壓翅片應該要特別注意驗算允許的壓降。

圖10 翅片厚度t對換熱性能的影響

3.2.4 翅片節距l的影響
圖11為h=7.1 mm、s=2.2 mm、t=6 mm時,不同l對j因子、f因子及FTEF因子的影響。由圖11可知,隨著l的增加,j因子、f因子、FTEF因子都降低,其中j因子降低30.1%,f因子降低41.5%,FTEF因子降低15.9%。翅片節距的增加,減小了對流體的擾動,j因子和f因子都減小。FTEF因子的減小,說明應該選擇短的翅片間距。

圖11 翅片節距l對換熱性能的影響

3.3 翅片參數的敏感性

圖12為翅片參數對于換熱器性能的敏感性。由圖可以發現,h和s對j因子呈正增長,而t和l對j因子呈負增長,且h和l對j因子的影響較大。h和t對f因子呈正增長,而s和l對f因子呈負增長。當t較大時,將會導致通道壓降快速增加,引起較高的運行費用。h和s對FTEF因子呈正增長,而t和l對FTEF因子呈負增長。這種從定性到定量分析翅片參數對換熱器性能影響的方法,對于指導翅片參數的選擇、提高板翅式換熱器性能具有重大的意義。

圖12 翅片參數對換熱器性能的敏感性

3.4 優化結果

圖13為以換熱器FTEF因子最大為目標函數的優化過程。從圖可知,在優化初始階段個體差異較大,較差的目標值被淘汰,目標函數值變化較明顯,隨后目標函數值變化幅度逐漸減弱,

最終趨于一個

穩定值。表1為維持入口Re=800時,以最大j因子、FTEF因子及最小f因子為目標的優化結果。由于基于響應面的優化方法,優化結果是由遺傳算法得到的近似解,所以為了驗證結果的準確性,對優化結果與通過CFD真實求解的結果進行了對比。由表1可知,優化值與CFD真實求解的值誤差很小,都在±10%以內,說明了利用Kriging響應面技術并應用遺傳算法可以有效地優化設計板翅式換熱器。在維持Re恒定時,得到最大j因子和最大FTEF因子的翅片結構參數非常相似,即采用薄的翅片厚度、高的翅片高度、短的翅片節距和適當大的翅片間距。為了減小翅片阻力,應該采用薄的翅片厚度、低的翅片高度、大的翅片間距和適當大的翅片節距。考慮到實際生產,在翅片通道Re=800時,得到鋸齒型板翅式換熱器性能最好的翅片參數為h=9.5 mm、s=2.2 mm、t=0.1 mm和l=3 mm。

圖13 以FTEF最大為目標函數的優化過程

項目h/mms/mmt/mml/mmjfFTEF最大j因子9.472.200.133.010.01920.05790.0496驗證值0.01900.06130.0482誤差/%0.9-5.52.8最小f因子4.702.880.137.750.01220.03500.0374驗證值0.01220.03820.0362誤差/%0.2-8.63.3最大FTEF因子9.472.170.113.010.01920.06580.0502驗證值0.01910.06090.0484誤差/%0.7-8.43.7

4 結 論

利用Kriging響應面技術,以最大j因子、FTEF因子和最小f因子為目標函數,應用遺傳算法優化設計了鋸齒型板翅式換熱器的翅片參數,可得出以下結論。

(1)鋸齒型翅片高度h與翅片間距s對換熱器綜合性能FTEF因子呈正增長,而翅片厚度t和翅片節距l呈負增長,其中翅片高度的影響為35%,翅片間距為45%,翅片厚度為-65%和翅片節距為-32%。

(2)在翅片通道Re=800時,鋸齒型板翅式換熱器的翅片高度h=9.5 mm、翅片間距s=2.2 mm、翅片厚度t=0.1 mm和翅片節距l=3 mm時性能最佳。

(3)基于Kriging響應面技術,應用遺傳算法可以有效地優化設計板翅式換熱器,并且該方法克服了傳統優化方法對于經驗關聯式的依賴,對于解決缺少經驗關聯式的復雜問題非常有效。

[1] SHAH R K, WEBB L. Compact and enhanced heat exchangers: heat exchangers theory and practice [M]. Washington, DC, USA: Hemisphere, 1983: 425-468.

[2] HOLLAND J H. Outline for a logical theory of adaptive systems [J]. Journal of the Association for Computing Machinery, 1962, 9(3): 297-314.

[3] PENG H, LING X. Optimal design approach for the plate-fin heat exchangers using neural networks cooperated with genetic algorithms [J]. Appl Therm Eng 2008, 28(5/6): 642-650.

[4] XIE G N, SUNDEN B, WANG Q W. Optimization of compact heat exchangers by a genetic algorithm [J]. Appl Therm Eng, 2008, 28(8/9): 895-906.

[5] 謝公南, 王秋旺. 遺傳算法在板翅式換熱器尺寸優化中的應用 [J]. 中國電機工程學報, 2006, 26(7): 53-57. XIE Gongnan, WANG Qiuwang. Geometrical optimization design of plate-fin heat exchanger using genetic algorithm [J]. Proceedings of the CSEE, 2006, 26(7): 53-57.

[6] 郭江峰, 許明田, 程林. 基于耗散數最小的板翅式換熱器優化設計 [J]. 工程熱物理學報, 2011, 32(5): 827-831. GUO Jiangfeng, XU Mingtian, CHENG Lin. Optimization design of plate-fin heat exchanger based on entransy dissipation number minimization [J]. Journal of Engineering Thermophysics, 2011, 32(5): 827-831.

[7] MISHRA M, DAS P K, SUNIL S. Second law based optimization of crossflow plate-fin heat exchanger design using genetic algorithm [J]. Appl Therm Eng, 2009, 29(14): 2983-2989.

[8] GHOLAP A K, KHAN J A. Design and multi-objective optimization of heat exchangers for refrigerators [J]. Applied Energy, 2007, 84(12): 1226-1239.

[9] HILBERT R, JANIGA G, BARON R, et al. Multi-objective shape optimization of a heat exchanger using parallel genetic algorithms [J]. Int J of Heat Mass

Transfer, 2006, 49(15/16): 2567-2577.

[10]NAJAFI H, NAJAFI B, HOSEINPOORI P. Energy and cost optimization of a plate and fin heat exchanger using genetic algorithm [J]. Appl Therm Eng, 2011, 31(10): 1839-1847.

[11]KAYS W M, LONDON A L. Compact heat exchanger [M]. 2nd ed. New York, USA: McGraw-Hill, 1964: 97-112.

[12]WIETING A R. Empirical correlations for heat transfer and flow friction characteristics of rectangular offset-fin plate-fin heat exchangers [J]. J Heat Transfer, 1975, 97(3): 488-490.

[13]JOSHI H M, WEBB R L. Heat transfer and friction in the offset-strip fin heat exchangers [J]. Int J of Heat Mass Transfer, 1987, 30(1): 69-84.

[14]MANSON S V. Correlations of heat transfer data and of friction data for interrupted plane fins staggered in successive rows, NACA tech note 2237 [R]. Washington, DC, USA: National Advisory Committee for Aeronautics, 1950.

[15]MOCHIZUKI S, YAGI Y, YANG W J. Transport phenomena in stacks of interrupted parallel-plate surfaces [J]. Exp Heat Transfer, 1987, 1(2): 127-140.

[16]DUBROVSKY E V, VASILIEV V Y. Enhancement of convective heat transfer in rectangular ducts of interrupted surfaces [J]. Int J of Heat Mass Transfer, 1988, 31(4): 807-818.

[17]MANGLIK R M, BERGLES A E. Heat transfer and pressure drop correlations for the rectangular offset strip fin compact heat exchanger [J]. Exp Therm Fluid Sci, 1995, 10(2): 171-180.

[18]KIM M S, LEE J, YOOK S J, et al. Correlations and optimization of a heat exchanger with offset-strip fins [J]. Int J of Heat Mass Transfer, 2011, 54(9/10): 2073-2079.

[19]YANG Y J, LI Y Z. General prediction of the thermal hydraulic performance for plate-fin heat exchanger with offset strip fins [J]. Int J of Heat Mass Transfer, 2014, 78(7): 860-870.

[20]WANG G G, SHAN S. Review of metamodeling techniques in support of engineering design optimization [J]. Journal of Mechanical Design, 2007, 129(4): 370-380.

[21]SACKS J, WELCH W J, MITCHELL T J, et al. Design and analysis of computer experiments [J]. Statistical Science, 1989, 4(4): 409-435.

(編輯 杜秀杰)

Optimization Design for Offset Fin in Plate Heat Exchanger with Genetic Algorithm

YANG Huizhu1,WEN Jian1,TONG Xin1,LI Ke1,LI Yanzhong1,WANG Simin2

(1. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China;2. School of Chemical Engineering and Technology, Xi’an Jiaotong University, Xi’an 710049, China)

In the current optimization for plate fin heat exchanger with offset fins, simultaneous multi-parameter optimization is rarely considered, or researches usually depend on empirical relations. A strategy by genetic algorithm (GA) combined with Kriging response surface is proposed, where the Kriging response surface provides an approximate relationship between the objective function and design variables. The fin heighth, fin spaces, fin thicknesstand interrupted lengthjof offset fins are taken as four optimization parameters, while maximumjfactor, minimumffactor and maximumFTEFfactor as three single objective functions at Reynolds number of 800. The results show that the effects of the fin height and fin space on the overall performanceFTEFfactor are positive and the effects of the fin thickness and the interrupted length are negative; the heat exchanger performance reaches the optimum as fin heighth=9.5 mm, fin spaces=2.1 mm, fin thicknesst=0.1 mm and interrupted lengthl=3 mm; the genetic algorithm combined with Kriging response surface eliminates the dependence on empirical relations.

plate-fin heat exchanger; genetic algorithm; Kriging response surface; optimization design

2015-04-27。

楊輝著(1988—),男,博士生;王斯民(通信作者),男,副教授。

國家自然科學基金資助項目(51106119,81100707);中央高校基本科研業務費專項資金資助項目。

時間:2015-10-28

10.7652/xjtuxb201512015

TK124

A

0253-987X(2015)12-0090-07

網絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20151028.1527.002.html

猜你喜歡
關聯優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
不懼于新,不困于形——一道函數“關聯”題的剖析與拓展
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
“一帶一路”遞進,關聯民生更緊
當代陜西(2019年15期)2019-09-02 01:52:00
奇趣搭配
智趣
讀者(2017年5期)2017-02-15 18:04:18
主站蜘蛛池模板: 国产精品专区第1页| 亚洲国产理论片在线播放| 婷婷亚洲最大| 久久先锋资源| 中文字幕日韩欧美| 一本久道热中字伊人| 亚洲日本www| 欧美另类第一页| 国产欧美日韩18| 精品超清无码视频在线观看| 人妻无码中文字幕一区二区三区| 成人噜噜噜视频在线观看| 日韩无码黄色| 国产成人在线无码免费视频| 国产剧情一区二区| 国产精品无码久久久久久| 亚洲一区二区成人| 在线永久免费观看的毛片| 精品夜恋影院亚洲欧洲| 免费毛片视频| 亚洲欧美国产高清va在线播放| 亚洲AV无码乱码在线观看裸奔| 精品人妻AV区| 无码国产伊人| 精品成人免费自拍视频| 久久国产免费观看| 综合色在线| 欧美成人二区| 2021国产精品自产拍在线| 亚洲无线国产观看| 色综合手机在线| 成年女人a毛片免费视频| 5555国产在线观看| 亚洲综合经典在线一区二区| 欧美三级不卡在线观看视频| 亚洲乱码在线播放| 91久久夜色精品国产网站| 亚洲区欧美区| 久久中文电影| 国产三级视频网站| 国产美女一级毛片| 性欧美精品xxxx| 国产全黄a一级毛片| 亚洲色图综合在线| 欧美国产视频| 国产91丝袜在线观看| 日本五区在线不卡精品| 九九九精品视频| 亚洲美女视频一区| 曰AV在线无码| 日本在线欧美在线| 国产女人18水真多毛片18精品| 免费毛片a| 青青国产视频| 丰满人妻久久中文字幕| 国产性生交xxxxx免费| 一本久道久久综合多人| 毛片大全免费观看| 欧美中文字幕在线视频| 伊人久久大香线蕉影院| 日韩最新中文字幕| 亚洲天堂免费观看| 日本在线国产| 国产呦视频免费视频在线观看| A级毛片无码久久精品免费| 亚洲欧洲一区二区三区| 国产精品美乳| 久久国产乱子| 国产成人三级| 美女啪啪无遮挡| 成人日韩视频| 黑人巨大精品欧美一区二区区| 狠狠色狠狠综合久久| 素人激情视频福利| 亚洲自偷自拍另类小说| 久久精品欧美一区二区| 欧美一级视频免费| 国产麻豆va精品视频| 欧美色香蕉| 亚洲高清在线播放| 欧美黄色网站在线看| 亚洲国产在一区二区三区|