亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        工程機械百葉窗翅片式散熱器的多目標(biāo)優(yōu)化設(shè)計

        2016-06-21 15:06:44張銀亮蔡惠坤
        中國工程機械學(xué)報 2016年6期
        關(guān)鍵詞:百葉窗翅片散熱器

        張銀亮, 蔡惠坤, 沈 超

        (1.廈門大學(xué) 機電工程系,福建 廈門 361005;2.廈門大學(xué) 深圳研究院,深圳,518057)

        工程機械百葉窗翅片式散熱器的多目標(biāo)優(yōu)化設(shè)計

        張銀亮1,2, 蔡惠坤1,2, 沈 超1

        (1.廈門大學(xué) 機電工程系,福建 廈門 361005;2.廈門大學(xué) 深圳研究院,深圳,518057)

        摘 要:為了提高工程機械百葉窗翅片式散熱器的散熱性能,提出了基于翅片質(zhì)量和綜合性能評價因子的多目標(biāo)優(yōu)化設(shè)計方法.在保證數(shù)值仿真方法可靠的前提下,在CFX中基于流固耦合的方法對翅片的傳熱特性和流阻特性進行了計算和分析.利用SolidWorks三維建模參數(shù)化模型后,在ANSYS Workbench仿真軟件中,運用多目標(biāo)遺傳算法對翅片質(zhì)量和綜合性能評價因子進行多目標(biāo)優(yōu)化,得到Pareto優(yōu)化解集.在保證翅片翅片綜合性能評價因子基本不變的前提下,通過優(yōu)化翅片質(zhì)量提高散熱器的散熱性能.結(jié)果表明:該方法適合百葉窗翅片的多參數(shù)優(yōu)化,在散熱器外形尺寸一定的情況下,提高了散熱器的散熱性能.綜合性能評價因子增加了1.94%,質(zhì)量減輕了21.00%.

        百葉窗翅片; 多目標(biāo)優(yōu)化; 散熱性能; 工程機械

        工程機械發(fā)動機的轉(zhuǎn)速和功率高,本身熱負(fù)荷很大,同時為了滿足工程機械NVH要求,使得噪聲等振動能量在封閉的發(fā)動機艙內(nèi)也被轉(zhuǎn)化成熱量.為了滿足工程機械的散熱性能需求,高效緊湊的百葉窗翅片式散熱器因其良好的傳熱特性和流阻特性,在車輛發(fā)動機散熱中得到廣泛應(yīng)用[1-2].

        為了進一步提高百葉窗翅片式散熱器的散熱效果,滿足在高原等惡劣環(huán)境中的工作要求,需要對百葉窗翅片進行結(jié)構(gòu)優(yōu)化.由于百葉窗翅片式結(jié)構(gòu)涉及多個參數(shù)(包括翅片間距、翅片厚度、翅片寬度、百葉窗角度、百葉窗間距),而利用三維軟件和仿真軟件的耦合,基于響應(yīng)面法的多目標(biāo)優(yōu)化方法在數(shù)控加工中心等結(jié)構(gòu)優(yōu)化中應(yīng)用成熟[3-6],適合對百葉窗翅片進行多參數(shù)的結(jié)構(gòu)優(yōu)化.

        本文針對某高原型挖掘裝載機冷卻液散熱器提出了一種以翅片間距、翅片厚度、翅片寬度、百葉窗角度、百葉窗間距為多參數(shù)輸入,以翅片質(zhì)量和綜合性能評價因子為多目標(biāo)的優(yōu)化設(shè)計方法.在CFX流體分析軟件中保證數(shù)值仿真方法可靠的前提下,對百葉窗翅片的傳熱特性和流阻特性進行了計算和分析后,在ANSYS Workbench仿真軟件中對百葉窗翅片進行了多目標(biāo)優(yōu)化設(shè)計.該方法適合百葉窗翅片的多參數(shù)結(jié)構(gòu)優(yōu)化,在散熱器外形尺寸一定的情況下,可以提高散熱器的散熱性能.

        1 計算模型的建立

        1.1 幾何模型參數(shù)

        表1和圖1某高原型挖掘裝載機四排百葉窗翅片式散熱器主要結(jié)構(gòu)尺寸.

        表1 百葉窗翅片尺寸Tab.1 Fin size

        圖1 計算模型Fig.1 Computational model

        1.2 控制方程及邊界條件

        使用ANSYS meshing劃分四面體網(wǎng)格,體網(wǎng)格最大尺寸為0.80,面網(wǎng)格最大尺寸為0.24.在CFX里面邊界條件的設(shè)定:

        (1) 空氣入口邊界為均勻速度空氣介質(zhì);

        (2) 空氣出口邊界設(shè)定為壓力出口條件;

        (3) 空氣與翅片接觸面設(shè)置成流固耦合傳熱面;

        (4) 水管壁面設(shè)置為恒溫,363 K;

        (5) 其他面設(shè)置為周期性對稱面.

        1.3 數(shù)據(jù)的處理

        工程機械百葉窗翅片式散熱器的傳熱特性和流阻特性主要取決于翅片結(jié)構(gòu)、管排數(shù)以及冷熱流體的熱物理性質(zhì)等參數(shù).其傳熱量和空氣側(cè)表面換熱系數(shù)基本的計算公式如下:

        圖2 網(wǎng)格的劃分Fig.2 Meshing

        (1)

        式中m——空氣質(zhì)量流量:kg/s;

        cpcp——水的定壓比熱:kJ/(kg·K);

        TinTin——空氣入口溫度:K ;

        Tout——空氣出口溫度:K.

        圖3 邊界條件Fig.3 Boundary condition

        百葉窗翅片空氣側(cè)對數(shù)平均溫差ΔT:

        (2)

        式中,Tw為水管壁面溫度(K).

        百葉窗翅片空氣側(cè)表面換熱系數(shù)h定義h如下

        (3)

        式中,Ao為空氣與翅片的換熱面積:(m2).

        在百葉窗翅片數(shù)值模擬計算過程中,空氣做不可壓縮假設(shè),其空氣側(cè)壓降可以直接通過數(shù)值模擬計算獲得.翅片的傳熱特性可以采用無量綱傳熱j因子來評價,其表達式為:

        (4)

        式中,Pr為空氣普朗特數(shù).

        翅片的流動阻力特性可以采用無量綱摩擦因子f來評價,其表達式為:

        (5)

        Δp=pin-pout

        (6)

        式中

        Δp——空氣側(cè)壓降(Pa);

        kc——冷卻空氣進出試件面積突縮產(chǎn)生的壓力損失系數(shù);

        ke——冷卻空氣進出試件面積突擴產(chǎn)生的壓力損失系數(shù);

        pin——空氣入口壓強(Pa);

        pout——空氣出口壓強(Pa).

        定義無量綱因子j/f1/3為散熱器綜合性能評價因子,對翅片性能做出綜合評價[4].

        1.4 計算模型的驗證

        為了判斷計算模型是否收斂,設(shè)置能量、質(zhì)量和速度殘差小于1×10-5時,最大迭代次數(shù)200步認(rèn)為結(jié)果收斂,同時監(jiān)控空氣出口溫度和壓降.如圖4所示,能量、質(zhì)量和速度殘差均小于1×10-5,迭代步數(shù)為80次.同時監(jiān)控的量空氣側(cè)壓降和出口溫度也達到了穩(wěn)定,說明計算結(jié)果已經(jīng)收斂.

        圖4 收斂及監(jiān)控量的平衡Fig.4 Convergence and balance of monitoring

        為了驗證計算模型的正確性,劃分五種不同網(wǎng)格數(shù)量,設(shè)置空氣入口風(fēng)速10 m/s,水管壁面溫度363 K,計算空氣進出口壓降和換熱系數(shù).如圖所示,當(dāng)網(wǎng)格數(shù)量達到129萬時,計算結(jié)果的偏差逐漸減小,在偏差小于3%~5%時,認(rèn)為網(wǎng)格數(shù)量對結(jié)果是可信的.綜合考慮,本計算模型選擇網(wǎng)格數(shù)量為253萬進行計算.

        為了驗證數(shù)值模擬的合理性,根據(jù)文獻的試驗條件,按照本文的數(shù)值模擬方法對百葉窗翅片散熱器進行了數(shù)值模擬,計算空氣進出口壓降和換熱系數(shù).數(shù)值模擬的結(jié)果和試驗數(shù)據(jù)能很好地吻合,空氣換熱系數(shù)在風(fēng)速較大時,偏差增大,最大偏差為12.23%,平均偏差為6.8%.數(shù)值模擬與試驗數(shù)據(jù)的偏差屬于可以接受的范圍之內(nèi),可以認(rèn)為計算模型合理,滿足工程的實際需要.

        圖5 網(wǎng)格獨立性驗證Fig.5 Grid independence verification

        圖6 計算結(jié)果的驗證Fig.6 Validation of computational results

        2 百葉窗翅片的多目標(biāo)優(yōu)化設(shè)計

        對于百葉窗翅片結(jié)構(gòu)的優(yōu)化,首先根據(jù)百葉窗翅片的結(jié)構(gòu)確定影響優(yōu)化目標(biāo)的關(guān)鍵參數(shù),在SolidWorks對百葉窗翅片建立三維模型;通過ANSYS Workbench(AWB)平臺與CAD接口技術(shù)連接,將參數(shù)化CAD模型轉(zhuǎn)化為AWB DS中的CAE模型;根據(jù)百葉窗翅片尺寸與干涉余量決策出參數(shù)值變化范圍.然后,采用中心復(fù)合試驗設(shè)計Central Composite Designs(CCD)確定試驗點,在AWB CFX仿真模塊中對試驗點進行流固耦合Fluid-Solid Coupling分析計算,根據(jù)多次試驗獲得的一組試驗數(shù)據(jù).

        在 ANSYS Workbench Design Explorer(AWB DE)中進行目標(biāo)驅(qū)動優(yōu)化Goal Driven Optimization(GDO)[8]模塊,也就是多目標(biāo)優(yōu)化技術(shù),從給出的一組樣本中得出“最佳”的設(shè)計點,同時觀察響應(yīng)曲線和響應(yīng)面的關(guān)系等.在n維可行解區(qū)域內(nèi)抽取均勻分布的樣本點,作為遺傳算法的初始種群[3-4].以綜合性能評價因子和質(zhì)量作為優(yōu)化目標(biāo),空氣側(cè)換熱系數(shù)和壓降作為約束,參數(shù)化結(jié)構(gòu)尺寸作為優(yōu)化設(shè)計變量,建立優(yōu)化模型.最后,采用多目標(biāo)優(yōu)化遺傳算法在決策空間中尋找到綜合性能評價因子大和翅片質(zhì)量小的優(yōu)化目標(biāo).

        2.1 響應(yīng)面模型的建立

        響應(yīng)面法[9](Response surface methodology,RSM)是一種試驗設(shè)計與數(shù)理統(tǒng)計相結(jié)合的優(yōu)化方法,其采用實驗設(shè)計理論對指定的設(shè)計點集合進行試驗,并在設(shè)計空間構(gòu)造測定量的全局逼近,得到目標(biāo)函數(shù)與約束函數(shù)的響應(yīng)面模型,來預(yù)測非試驗點的響應(yīng)值.

        為了建立百葉窗翅片的響應(yīng)面模型,必須先確定翅片的設(shè)計變量.由于百葉窗翅片結(jié)構(gòu)復(fù)雜,其散熱性能受結(jié)構(gòu)尺寸影響,通過理論分析和查閱相關(guān)資料,將對其影響較明顯的五個結(jié)構(gòu)尺寸作為設(shè)計變量:翅片厚度,翅片間距,翅片寬度,百葉窗間距,百葉窗角度.

        依據(jù)試驗設(shè)計理論和經(jīng)驗,確定百葉窗翅片選取的5個設(shè)計變量的取值范圍,通過中心復(fù)合試驗設(shè)計[10]在其組成的決策空間中構(gòu)造出27個試驗設(shè)計點.

        2.2 數(shù)學(xué)模型的建立

        在工程機械散熱器的多目標(biāo)優(yōu)化設(shè)計中,一般都是在散熱器外形尺寸一定的情況下,以空氣側(cè)壓降和換熱系數(shù)為目標(biāo)函數(shù).而在數(shù)值仿真過程中,由于工程機械散熱器外形尺寸較大,只能以部分翅片單元進行數(shù)值仿真.所以在本文散熱器優(yōu)化設(shè)計過程中,以翅片綜合性能評價因子和翅片質(zhì)量為目標(biāo)函數(shù).而對于百葉窗翅片結(jié)構(gòu)的優(yōu)化問題,總是希望翅片綜合性能評價因子大,同時翅片質(zhì)量小,從而在散熱器外形尺寸一定的情況下,提高散熱器的散熱性能.

        基于上述建立的響應(yīng)面模型,以翅片綜合性能評價因子和質(zhì)量為多目標(biāo)優(yōu)化的目標(biāo)函數(shù),翅片空氣側(cè)壓降和換熱系數(shù)作為約束條件,翅片的五個結(jié)構(gòu)尺寸為設(shè)計變量.得到翅片目標(biāo)優(yōu)化數(shù)學(xué)模型如下:

        min(-jf);min(m);s.t.h≥;

        (i=1,2,…,n)

        (7)

        式中:

        X——決策矢量,X∈Ω為決策空間;

        di——設(shè)計尺寸變量;

        jf——綜合性能評價因子;

        h——空氣側(cè)換熱系數(shù);

        p——空氣側(cè)壓降;

        2.3 Pareto最優(yōu)解的獲得

        在單目標(biāo)優(yōu)化中,往往能獲得一個最優(yōu)解,但在多目標(biāo)優(yōu)化中,各目標(biāo)之間很難同時達到最優(yōu),所以多目標(biāo)優(yōu)化常常產(chǎn)生一系列Pareto解,也叫做有效解.

        求解百葉窗翅片的多目標(biāo)優(yōu)化問題就是盡可能找到更多具有代表性的符合約束條件的Pareto解,在計算獲得平均分布的Pareto解之后,根據(jù)設(shè)計要求和工程經(jīng)驗,選擇最合理的百葉窗翅片結(jié)構(gòu)參數(shù).本文選用多目標(biāo)遺傳算法對式(7)的數(shù)學(xué)模型進行優(yōu)化來得到百葉窗翅片的全局Pareto解.

        2.4 優(yōu)化結(jié)果

        根據(jù)原有工程機械散熱器翅片尺寸的初步設(shè)計,將CAD模型導(dǎo)入AWB中,對五個核心尺寸進行參數(shù)化定義,并進行網(wǎng)格劃分.然后設(shè)置邊界條件,初步驗證網(wǎng)格獨立性和仿真模型的準(zhǔn)確性,初步仿真得到綜合性能評價因子為0.013 276,質(zhì)量為0.002 104 1 kg,空氣側(cè)壓降為997.93Pa,換熱系數(shù)為96.51 W/(m2·K).最后在此基礎(chǔ)上利用AWB DE目標(biāo)驅(qū)動優(yōu)化進行結(jié)構(gòu)優(yōu)化,得到多目標(biāo)的Pareto最優(yōu)解,如圖7所示,選取其中六組解如表2所示.

        表2 多目標(biāo)優(yōu)化解集Tab.2 Multi-objective optimal solution set

        圖7 多目標(biāo)的Pareto圖Fig.7 Multi objective Pareto graph

        為了得到最終優(yōu)化的翅片結(jié)構(gòu),需要對參數(shù)進行靈敏度分析.如圖8所示,五個設(shè)計參數(shù),從左到右依次是翅片間距、翅片厚度、翅片寬度、百葉窗角度、百葉窗間距.從各結(jié)構(gòu)參數(shù)的靈敏度可以看出,翅片間距和翅片寬度與綜合性能評價因子成負(fù)相關(guān),翅片厚度和百葉窗間距與綜合性能評價因子成正相關(guān),百葉窗角度對綜合性能評價因子影響很小.翅片間距、翅片厚度、翅片寬度與質(zhì)量成正相關(guān),而百葉窗間距與綜合性能評價因子成正相關(guān),百葉窗角度和百葉窗間距對質(zhì)量的影響很小.因此在這里需要減小翅片間距和翅片厚度,增大百葉窗間距.

        圖8 靈敏度分析Fig.8 Sensitivity analysis

        如下表3所示給出了原始設(shè)計參數(shù)和優(yōu)化后的設(shè)計參數(shù)與優(yōu)化前后結(jié)果對比.從原始設(shè)計翅片性能與最終優(yōu)化取整后翅片性能比較后可知,壓降增加了1.91%,換熱系數(shù)減少了1.41%,綜合性能評價因子增加了0.60%,質(zhì)量減輕了20.08%.取整后的參數(shù),可看出換熱系數(shù)變化略大減少了4.1%,但是壓降減少了2.29%,綜合性能評價因子增加了1.94%,質(zhì)量減輕了21.00%.從綜合性能評價因子的角度來看的話,翅片整體的散熱性能是沒有太大變化,但是質(zhì)量卻減輕了21%,這樣在散熱器外形尺寸一定的情況下,通過增加散熱器的散熱面積,從而能夠提高散熱器的散熱性能.

        3 結(jié)論

        (1) 建立了工程機械百葉窗翅片的3D參數(shù)化模型,基于AWB DE 仿真優(yōu)化模塊轉(zhuǎn)化為有限元模型.利用CFX流體數(shù)值仿真分析,并對數(shù)值仿真模型進行了驗證.通過分析百葉窗翅片的散熱特性和流阻特性,為百葉窗翅片的多目標(biāo)優(yōu)化提供參考模型.

        表3 優(yōu)化參數(shù)前后結(jié)果對比Tab.3 Comparison of the results before and after optimization

        (2) 基于響應(yīng)面模型,得到了工程機械百葉窗翅片的綜合性能評價因子和質(zhì)量的靈敏度分析圖,以及百葉窗翅片全局Pareto最優(yōu)解.

        (3) 本文將流固耦合分析、響應(yīng)面模型、抽樣技術(shù)、多目標(biāo)遺傳算法和靈敏度分析法結(jié)合,對工程機械百葉窗翅片的結(jié)構(gòu)參數(shù)進行了優(yōu)化.最終結(jié)果表明:綜合性能評價因子增加了1.94%,質(zhì)量減輕了21.00%.

        [1] 漆波,李隆鍵,崔文智,陳清華.百葉窗式翅片換熱器中的耦合傳熱[J].重慶大學(xué)學(xué)報,2005,28(10):39-42.

        QI Bo,LI Longjian,CUI Wenzhi,CHEN Qinghua.Coupled conduction convective heat transfer in the louvered fin heat exchanger[J].Journal of Chongqing University,2005,28(10):39-42.

        [2] 董軍啟,陳江平,陳芝久.百葉窗翅片的傳熱與阻力性能試驗關(guān)聯(lián)式[J].制冷學(xué)報,2007,28(5):10-14.

        DONG Junqi,CHEN Jiangping,CHEN Zhijiu.Correlation of flow and heat transfer characteirstics for multi-louvered fin[J].Junrnal of Refrigeration,2001,28(5):10-14.

        [3] 姜衡,管貽生,邱志成.基于響應(yīng)面法的立式加工中心動靜態(tài)多目標(biāo)優(yōu)化[J].機械工程學(xué)報,2011,47(11):125-133.

        JIANG Heng,GUAN Yisheng,QIU Zhicheng.Dynamic and static multi-objective optimization of a vertical machining center based on responsesurface method[J].Journal of Mechanical Engineering,2011,47(11):125-133.

        [4] 錢堯一,侯亮,卜祥建.發(fā)動機艙排氣引射系統(tǒng)多目標(biāo)優(yōu)化設(shè)計[J].工程機械,2014,45(1):40-45.

        QIAN Yaoyi,HOU Liang,BU Xiangjian.Multi-objective optimal designof underhood exhaust ejector system in engine compartment[J].Construction Machinery and Equipment,2014,45(1):40-45.)

        [5] 覃祖和,王志越,黃美發(fā),林振廣.數(shù)控銑床工作臺多目標(biāo)優(yōu)化研究[J].機械設(shè)計與制造,2016(1):101-104.

        QIN Zuhe,WANG Zhiyue,HUANG Meifa,LIN Zhenguang.Research on multi-objective optimization of numerical control milling machine worktable[J].Machinery Design & Manufacture,2016(1):101-104.

        [6] 毛航,王珂,王永慶,劉彤,劉敏珊.基于響應(yīng)面分析法的百葉窗翅片結(jié)構(gòu)優(yōu)化設(shè)計[J].化工設(shè)備與管道,2015,52(1):23-27.

        MAO Hang,WANG Ke,WANG Yongqing,LIU Tong,LIU Minshan.Optimum design of louvered fin structure based on response surface methodology[J].Process Equipment & Piping,2015,52(1):23-27.

        [7] 朱家玲,李曉光,張偉.車用散熱器百葉窗布置方式的數(shù)值模擬與分析[J].天津大學(xué)學(xué)報,2013,46(3):244-249.

        ZHU Jialing,LI Xiaoguang,ZHANG Wei.Numerical simulation and analysis of different louver arrangements of automotive heat exchangers[J].Journal of Tianjin University,2013,46(3):244-249.

        [8] 蒲廣益.基于ANSYS Workbench12基礎(chǔ)教程與實例詳解[M].中國水利水電出版社,2010:244.

        PU Guangyi.Tutorials and detailed examples of ANSYS Workbench12[M].China Water & Power Press,2010:244.

        [9] Myers R H,Montgomery D C.Response surface methodology:Process and product optimization using designed experiments[M].Wiley,1995:1-78,351-401.

        [10] BOX G E P,DRAPER N R.Experimental model buildingand response surfaces[M].John Wiley & Sons,1987.

        Multi-objective Optimization for Engineering Machinery Radiator with Louvered Fins

        ZHANG Yin-liang1,2, CAI Hui-kun1,2, SHEN Chao1

        (1.School of Areospace Engineering,Xiamen University,Xiamen 361005,China; 2.Shen Zhen Research Institute,Xiamen University,Shenzhen 518057,China)

        In order to improve the heat transfer performance of construction machinery radiator with louvered fins, a multi-objective optimization design method is proposed in terms of mass and comprehensive evaluation factors. With the reliability of numerical simulation method, the fin heat transfer and flow resistance characteristics are calculated via CFXTM based on fluid-structure coupling. By establishing three-dimensional parametric model using SolidWorksTM and ANSYS WorkbenchTM, the higher comprehensive evaluation factor and lighter weight are obtained, whereas the Pareto optimal solution set is attained by the multi-objective genetic algorithm for optimization. Under the same comprehensive performance evaluation factor of fin, the fin heat dissipation performance is improved by optimizing the fin quality. Therein, it is detected from results that the approach is suitable for multi-parameter optimization on louvered fins, whilst the radiator heat transfer performance is improved using the fixed size of radiator. In addition, the comprehensive performance evaluation factor is increased by 1.94%, whereas the quality is reduced by 21%.

        louvered fin; multi-objective optimization; heat transfer performance; construction machinery

        福建省自然科學(xué)基金項目(2014J01210);深圳市科技計劃項目(JCYJ2014041716249675)

        張銀亮(1991-),男,碩士.E-mail:Zhang_yin_liang@163.com

        TH 123,U 464.138.2

        A

        1672-5581(2016)06-0508-07

        猜你喜歡
        百葉窗翅片散熱器
        垂直翅片管自然對流傳熱特性的數(shù)值研究
        機械工程師(2022年6期)2022-06-21 08:44:24
        ◆ 散熱器
        散熱器
        大功率COB-LED的翅片散熱器優(yōu)化研究
        ◆ 散熱器
        散熱器
        超硬翅片滾刀加工在CNC磨床上的實現(xiàn)
        讓百葉窗動起來FlipFlic百葉窗開啟器
        發(fā)電百葉窗
        縱向渦發(fā)生器對百葉窗翅片管換熱器性能的提升
        亚洲av理论在线电影网| 综合图区亚洲另类偷窥| 亚洲国产精品国自产拍av| 欧美白人最猛性xxxxx| 亚洲最稳定资源在线观看| 国产亚洲精品一区二区在线观看| 97在线视频人妻无码| 一本无码人妻在中文字幕免费| 欧洲日韩视频二区在线| 亚洲国产综合精品一区最新| 亚洲av无码乱码国产麻豆| 999久久久国产精品| 香蕉国产人午夜视频在线观看| 国产精品自拍网站在线| 免费不卡无码av在线观看| 日韩a无v码在线播放| 久久中文字幕亚洲精品最新| 少妇被猛烈进入中文字幕| 美女不带套日出白浆免费视频| 亚洲成人小说| 青青草免费激情自拍视频| 日本久久久免费观看视频| 欧美寡妇xxxx黑人猛交| 蜜桃av中文字幕在线观看| 色综合中文综合网| 国产av精国产传媒| 黑人性受xxxx黑人xyx性爽| 毛片在线啊啊| 国产一区二区三区不卡在线播放| 日本一二三区在线观看视频| 99国产精品自在自在久久| 午夜无码片在线观看影院| 国产精品福利久久香蕉中文| 久久国产精品美女厕所尿尿av| 久久精品亚洲精品国产色婷| 特黄做受又粗又长又大又硬| 精品国产a∨无码一区二区三区| 丝袜 亚洲 另类 欧美| 亚洲国产精品国自拍av| 日韩人妻一区二区三区蜜桃视频| 亚洲免费不卡|