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

        ?

        基于設計空間探索的型線自動優(yōu)化

        2010-06-07 07:53:46許勁松楊小玉
        船舶力學 2010年7期
        關鍵詞:型線潛艇阻力

        梁 軍,許勁松,謝 杰,楊小玉

        (1上海交通大學船舶海洋與建筑工程學院,海洋工程國家重點實驗室,上海 200030;2上海馮卡門計算機科技有限公司,上海 200235)

        基于設計空間探索的型線自動優(yōu)化

        梁 軍1,許勁松1,謝 杰1,楊小玉2

        (1上海交通大學船舶海洋與建筑工程學院,海洋工程國家重點實驗室,上海 200030;2上海馮卡門計算機科技有限公司,上海 200235)

        船舶型線優(yōu)化是改善船舶阻力性能、提高運行經濟性的關鍵環(huán)節(jié)。文中將型線修改、阻力CFD模擬、設計空間探索優(yōu)化等流程整合在一起,創(chuàng)建了一個型線“自動優(yōu)化”的平臺。在該平臺上實現(xiàn)了對美國海軍水面作戰(zhàn)中心(Naval Surface Warfare Center)潛艇模型DTRC Model 5470艏部型線的自動優(yōu)化,得出了一個阻力性能最優(yōu)的艏部型線方案。經過模型拖曳試驗驗證,證明了該種“自動優(yōu)化”方案的可行性與有效性。

        自動優(yōu)化;船舶型線;阻力性能

        1 引 言

        船舶阻力性能是影響船舶運行經濟性的重要因素,通過對船舶線型的優(yōu)化可以大幅提高船舶的節(jié)能性,對世界海運業(yè)的健康發(fā)展具有舉足輕重的作用[1]。

        傳統(tǒng)的型線優(yōu)化設計大多按照經驗修改、CFD模擬計算、性能評價的步驟重復計算多個方案,然后在這些方案中進行選擇。在美國一項三體船阻力性能優(yōu)化的課題研究中[2],對外船體與中心船體的縱向位置和橫向間距的多種組合進行了消波效果探討,共完成了48 678個組合方案的CFD計算,從中選擇確定了最佳設計方案。這種優(yōu)化過程可稱為“手工優(yōu)化”[3]。

        “手工優(yōu)化”方法雖然已大量應用,但無法突破其內在的局限性。第一、需要多次手動修改原始線型方案,反復手工設置基本一致的CFD模擬,然后對計算所得結果進行比較。重復性的人工干預過程耗費了大量的資源,也造成了計算機的間斷性工作,無法大幅提高優(yōu)化工作效率。第二、“手工優(yōu)化”所獲得的優(yōu)化方案局限于設計者的經驗,常常只能在一定程度上改進阻力性能,而不能保證獲得全設計空間中的最優(yōu)值。這種優(yōu)化設計本質上只是方案改進,而不是真正意義上的最優(yōu)設計。

        為克服“手工優(yōu)化”的不足,人們一直在探索將CFD模擬與自動優(yōu)化方法相結合的“自動優(yōu)化”平臺[3-6],期望達成下列目標:第一、充分利用計算機做重復有序性工作的優(yōu)越能力,在多種軟件信息交流的基礎上,將整個設計優(yōu)化過程集成為一個自動流程,從而達到節(jié)省資源、提高效率的目標;第二、基于最優(yōu)理論,綜合考慮各約束條件,對設計空間進行深入而全面的探索,以獲得全設計空間中的最優(yōu)方案。

        本研究將型線修改、阻力CFD模擬、設計空間探索優(yōu)化等流程整合在一起,利用OPTIMUS 5.2將GAMBIT2.3.16及FLUENT6.3.26中的計算過程集成起來,創(chuàng)建了一個“自動優(yōu)化”的平臺。在該平臺上實現(xiàn)了對美國海軍水面作戰(zhàn)中心(Naval Surface Warfare Center)潛艇模型DTRC Model 5470[7-8]艏部型線的自動優(yōu)化,得出了一個阻力性能最優(yōu)的艏部型線方案。經過模型拖曳試驗驗證,證明了該種“自動優(yōu)化”方案的可行性與有效性。

        2 型線“自動優(yōu)化”方法

        2.1 基于FLUENT的潛艇阻力CFD計算

        本研究以潛艇模型DTRC Model 5470為案例[7-8],計算模型縮尺比為1:50,模型總長2.178m,計算模型外形見圖1。

        計算采用無量綱化的RANS方程作為控制方程,張量形式表達為:

        圖1 潛艇模型DTRC Model 5470Fig.1 Submarine model DTRC Model 5470

        其中,ρ為流體的密度,υ為流體的運動粘性系數,F(xiàn)i為外力項,ui為平均速度,p為平均壓力,ui′為速度脈動量。雷諾應力項借助湍流模型表達以實現(xiàn)方程的封閉:

        其中,δij是Kronecker符號,vt為紊動粘性系數,Cμ是常數;k為紊動能及其耗散率為ε。

        湍流模式采用RNG k-ε兩方程模型,通過在大尺度運動和修正后的粘度體現(xiàn)小尺度的影響,使小尺度運動有系統(tǒng)地從控制方程中去除,可以更好地處理高應變率及流線彎曲程度較大的流動模擬,

        其中,Gk是由層流速度梯度而產生的湍流動能,Gb是由浮力而產生的湍流動能,YM是在可壓縮湍流中因過渡的擴散產生的波動,C1ε=1.42,C2ε=1.68是常量,αk和αε分別是k方程和ε方程的湍流Prandtl數,Sk和Sε可由用戶定義為0。

        在GAMBIT2.3.16中生成三維非結構化網格。相對于結構化網格,非結構化網格的適應能力強,局部加密也比較容易。在梯度大的地方,網格必須保證足夠細密。為了有效而接近實際的壁面模擬,要求合理布置邊界層區(qū)域內的網格[9]。定義量綱一的量y+=uτΔyp/ν,其中Δyp為第一層網格結點離開壁面的距離,可由下式進行估算:

        較好的網格劃分應使y+在30到500之間,而且邊界層內有足夠的網格數。本案例計算中,Re=4.8×106(以潛艇長為特征長度),y+≈33.05。計算流場區(qū)域設定為高7m,半徑為2m的半圓柱體,潛艇首部離來流入口約一個艇長。為了減少網格數目并降低計算量,在潛艇外劃分一個長3m、半徑0.5m的小半圓柱區(qū)域,使網格總數控制在380萬左右。計算區(qū)域對稱面網格如圖 2所示。

        本案例計算采用定常計算方法,邊界條件設定如下:

        ●速度入口:潛艇艏部向前2m,x方向來流速度大小為航速,y和z方向為零;

        ●壁面:潛艇外表面;

        ●壓力出口:潛艇艉部向后約兩個艇長,設定相對參考壓力點的流體靜壓值;

        ●對稱面:垂直于對稱面的速度分量為0,平行于對稱面的速度分量的法向導數為0;

        ●外場:距潛艇表面約為2m,速度為未受到擾動的主流區(qū)速度。

        2.2 基于OPTIMUS的優(yōu)化方法

        由于目前應用較多的“手工優(yōu)化”方法主要根據模擬計算結果進行經驗修正,所獲得的優(yōu)化結果能夠優(yōu)于原型方案但無法保證最優(yōu)。最優(yōu)解的獲得必須建立在全設計空間探索的基礎上,一般包括試驗設計、響應面模擬、最優(yōu)點搜索等若干步驟。

        試驗設計DOE(Design Of Experiment)的目的是對數值或物理試驗進行科學合理的安排,以較少的試驗次數獲得較多的設計空間信息,達到設計空間探索的最佳效果。一個科學而完善的試驗設計,能夠合理地安排各種試驗因素,有效地分析試驗數據,從而使用較少的資源最大限度地獲得豐富而可靠的資料。試驗設計學科由費希爾在農業(yè)試驗研究中創(chuàng)立,運用均方差排列的拉丁方和方差分析原理解決了長期存在的實驗條件不均衡問題。二戰(zhàn)以后試驗設計已經成為不同領域各類試驗的通用技術[10]。

        在OPTIMUS中包含有近二十種試驗設計方法,如全因子(Full Factorial)、拉丁超立方體(Latin Hypercube)、隨機法(Random)、部分因子(Fractional Factorial)等。本次設計選用拉丁超立方體作為DOE方法,是一種“充滿空間”(space filling)的試驗設計方法。依據此方法進行若干數值模擬試驗,可獲得設計空間中數目與位置確定的一系列設計點,在此基礎上可以進行目標變量與控制變量的相關性分析等后處理,并且能夠以回歸方法建立控制變量與目標變量之間的函數關系,即響應面模型。

        響應面模型RSM(Response Surface Model)反映了目標變量(因變量)與多個控制變量(自變量)間的函數關系。由于這種函數關系一般是曲線或曲面的關系,因而稱為響應面模型RSM。由于響應面模型必須根據系列試驗數據回歸得到,回歸分析的優(yōu)劣程度直接決定了響應面模型的精確性?;貧w分析的過程一般可分為兩個階段,第一階段的主要目的是確定當前的設計點或試驗點是否接近響應面的最優(yōu)(最大或最?。┪恢?。當試驗點遠離響應面的最優(yōu)位置時,可使用如下的一階模型(first-order model)去逼近:

        圖2 非結構化網格Fig.2 Non-Structure mesh

        其中βi表示xi的斜率或線性效應。當試驗區(qū)域接近響應面的最優(yōu)區(qū)域或者位于最優(yōu)區(qū)域中時,可以開始第二階段的設計,目的就是獲得響應面在最優(yōu)值周圍的一個精確逼近并且識別出最優(yōu)設計點。這時常采用如下的二階模型(Second-order model)來逼近[11]:

        OPTIMUS為RSM響應面模擬提供了一系列的數學近似方法,如線性插值(Interpolation)、泰勒法(Taylor)等。在合適的數據量基礎上,使用Taylor法能得到近似程度相當不錯的響應面模型。而數據比較少時,可采用線性插值建立響應面模型。本次案例計算采用線性插值方法中的Kriging函數法[11]。

        在生成的RSM響應面模型上進行最優(yōu)點搜索,通常需要包括以下設定:設計目標(最大或最?。⒓s束條件、以及最優(yōu)算法的選擇。通過OPTIMUS能夠方便地給模型添加數學約束、確定設計目標,并提供了一系列的最優(yōu)點搜索算法,如梯度算法、序列二次規(guī)劃算法(Sequential Quadratic)等局部搜索算法,以及自適應遺傳算法(Self-Adaptive Evolution)、模擬退火算法等全局搜索算法。本次案例優(yōu)化采用自適應遺傳算法,是演化算法的一個分支[13]。

        2.3 “自動優(yōu)化”過程集成

        過程集成是將某個系統(tǒng)中各個獨立的流程單元通過良好的數據交換接口彼此連接,并按照整個系統(tǒng)的運行順序整合起來,形成一個完整的流程集成平臺。隨著計算機應用技術的提高,過程集成方法逐步滲透到了多種領域,尤其是設計制造領域,可以實現(xiàn)產品設計生產的高度自動化,達到節(jié)省資源、提高效率的目的。如波音公司將CAD建模、CFD性能模擬等過程集成,實現(xiàn)了設計、分析、加工和檢測的一體化,成功完成了波音777客機的設計生產,成為過程集成方法的應用典范。

        本次案例研究將潛艇型線修改、阻力CFD模擬、設計空間探索優(yōu)化等流程整合在一起,利用OPTIMUS 5.2將GAMBIT2.3.16及FLUENT6.3.26中的計算過程集成起來,創(chuàng)建了一個“自動優(yōu)化”的平臺。具體流程控制過程為:在GAMBIT中生成潛艇幾何模型,并定義型線修改的控制變量;在GAMBIT中生成計算模型并進行網格劃分,導出mesh文件至FLUENT;在FLUENT中進行模擬計算,獲得阻力計算結果;在OPTIMUS中設置優(yōu)化流程,根據試驗設計DOE結果改變型線控制變量,輸入GAMBIT;在GAMBIT中按照新的控制變量重新生成模型,重復上述流程,直至獲取DOE要求的足夠的數值試驗結果;在OPTIMUS中建立RSM響應面模型,完成最優(yōu)方案搜索。完整的流程結構如圖3所示。

        為了在OPTIMUS中實現(xiàn)上述流程,必須建立如圖4所示的過程集成平臺,包含下列流程單元:

        ●控制變量:通過該變量的數值變化實現(xiàn)潛艇型線和幾何模型的變化??稍贕AMBIT命令記錄文件gambit.jou中進行變量設置,達到自動修改的目的。在本次案例計算中,選取了潛艇子午線首部三個點的位置偏移量作為控制變量,可以控制首部型線修改;

        圖3 流程結構圖Fig.3 Flow of the procedures

        ●輸入文件gambit.jou,fluent.jou:gambit.jou為GAMBIT命令記錄文件,按照該文件運行可使GAM-BIT自動重復幾何模型生成、網格劃分、網格文件輸出的過程。fluent.jou為FLUENT命令記錄文件,按照該文件運行可使FLUENT自動重復讀入網格文件、邊界條件設定、完成模擬計算、輸出計算結果、保存case和data文件的過程;

        ●求解器FLUENT,GAMBIT:兩個求解器代表的是安裝在Linux服務器上的GAMBIT和FLUENT程序。OPTIMUS中通過添加遠程調用命令(rsh)來實現(xiàn)遠端服務器上的兩個程序的調用運行。GAMBIT以gambit.jou為輸入,輸出網格文件marin.msh;FLUENT以fluent.jou為輸入,讀入網格文件marin.msh后輸出結果文件out.trn。兩個求解器都包含有OPTIMUS自帶的等待命令,以保證輸出文件的正常生成和流程的連貫執(zhí)行。另外,為不使自動迭代過程中不同數值試驗生成的多個msh、cas、dat文件發(fā)生沖突,使用dos下的del命令將每一次執(zhí)行所生成的結果文件刪除,同時OPTIMUS將自動進行讀取和保存數值試驗結果。

        ●輸出文件out.trn:out.trn中包含F(xiàn)LUENT計算結果。由于out.trn格式固定,OPTIMUS可通過位置匹配將目標變量值和約束變量值抽取出來。

        ●輸出變量:OPTIMUS按照固定抽取規(guī)則從out.trn中獲得輸出變量,包括目標變量和約束變量。本次案例計算中抽取潛艇阻力(resistance)作為目標變量。

        圖4 OPTIMUS過程集成流程圖Fig.4 Integrated procedures in OPTIMUS

        在過程集成完成后,將控制變量的值設為零,完成一次試運行,測試流程的正確性以及模擬計算的準確性,這個過程在OPTIMUS中稱為Nominal,是一個流程調試過程。在本案例優(yōu)化中,以DTRC Model 5470潛艇原型的模擬計算作為流程調試Nominal依據,通過模型拖曳試驗進行驗證。

        3 案例優(yōu)化結果

        3.1 案例原型模擬結果與驗證

        DTRC Model 5470潛艇原型的模擬計算在小型計算服務器上完成,一個速度計算需要完成迭代過程1 000次,6核CPU并行計算約4小時可以收斂。

        原型模擬共計算了4個速度的阻力值,同時在上海交通大學船模拖曳水池完成了同尺度模型的拖曳試驗(見圖 5)。模擬計算與拖曳試驗的對比結果列于表1,數據吻合度較高,模擬算法的準確性能夠滿足“自動優(yōu)化”的要求。

        圖5 潛艇原型拖曳試驗Fig.5 Model test of original hull

        表1 原型潛艇阻力計算值與試驗值比較Tab.1 Comparison of computational and experimental resistance results

        3.2 自動優(yōu)化結果

        在本次案例優(yōu)化中設定了如圖6所示的3個控制變量:第一個控制變量offset_1為潛艇前端頂點x=0在x方向的偏移量,變動范圍為-0.1~0.01m;第二個控制變量offset_2為子午線上點x=0.035m沿y方向的偏移量,變動范圍為±0.02~±0.015m;第三個控制變量offset_3為子午線上點x=0.115m沿y方向的偏移量,變動范圍為±0.02~±0.015m。

        圖6 控制變量Fig.6 Control variables

        控制變量的組合方案由試驗設計DOE確定。本次優(yōu)化過程采用了可選擇試驗次數的拉丁方方法進行試驗設計,設定試驗次數為16次。OPTIMUS以隨機種子值為基礎,在輸入變量的區(qū)間范圍內以拉丁方算法的規(guī)則尋找16個組合方案,充滿設計空間(space filling),并按照這16個組合方案自動重復執(zhí)行4m/s速度下的數值試驗流程,獲得16個數值試驗結果,列于表2中。本案例DOE過程在6核CPU并行環(huán)境下連續(xù)自動運行約70小時后完成。

        表2 DOE16次模擬試驗結果Tab.2 16 DOE simulation results

        在DOE所獲得的數值試驗數據基礎上,可回歸建立響應面模型(RSM),獲得整個設計空間中控制變量與目標變量之間的函數關系模型。本次優(yōu)化過程采用Kriging插值法來建立響應面模型,其中的一個RSM面如圖7所示。

        本次優(yōu)化是阻力最小的單目標優(yōu)化,使用自適應進化算法在RSM面上搜尋最優(yōu)方案,獲得速度4m/s時的最優(yōu)阻力值為50.28N。該最優(yōu)方案對應的三個控制變量offset_1,offset_2,offset_3的值分別為-0.080 1m,0.011 3m,0.007 8m,優(yōu)化前后的潛艇型線對比如圖8所示,實物模型對比見圖9。

        3.3 優(yōu)化結果驗證

        為了驗證優(yōu)化結果的準確性,首先進行了CFD模擬驗證,根據最優(yōu)方案對應的控制變量結果重新生成計算模型,模擬計算結果為總阻力50.37N,與RSM上的優(yōu)化搜尋結果50.28N吻合較好。

        隨后在上海交通大學船模拖曳水池進行了優(yōu)化模型的拖曳試驗(見圖 10),試驗結果與計算結果列于表 3中,在3m/s、4m/s、5m/s都獲得了較高的一致性,只有2m/s時誤差大于10%。

        圖7 控制變量1、2與阻力的RSM面Fig.7 RSM of control variables 1,2

        表3 優(yōu)化潛艇阻力計算值與試驗值比較Tab.3 Comparison of computational and experimental resistance results

        將優(yōu)化前后的總阻力系數計算值和試驗值在圖11中同時繪出無因次曲線,可以看出3m/s至4m/s的中速區(qū)阻力減小2%左右,而在2m/s和5m/s附近區(qū)域阻力降低不多。這種優(yōu)化效果的差別是由于優(yōu)化過程針對4m/s速度點,對低速和高速區(qū)域不具有針對性。

        4 結 論

        本次研究通過對案例潛艇的優(yōu)化應用和試驗驗證,證明了所構建的型線“自動優(yōu)化”平臺方案的可行性。優(yōu)化過程中通過DOE獲得的16個試驗方案較好地反映了全設計空間的完整信息,并在70小時內實現(xiàn)了無人工干預的型線自動優(yōu)化,達到了節(jié)省資源、提高效率的目標。

        從計算和試驗結果來看,阻力變化趨勢十分吻合,優(yōu)化效果也很明顯。由于案例研究選用的DTRC Model 5470潛艇本身具有很好的阻力性能,因此阻力再優(yōu)化難度很大,優(yōu)化結果只取得了2%左右的阻力改進。

        下一步研究將把該“自動優(yōu)化”平臺進一步完善后應用于水面船舶的阻力性能優(yōu)化,尋求在兩個方面克服瓶頸問題:第一、解決好船體NURBS曲面輸出轉化問題,將船舶型線設計的CAD工具集成入“自動優(yōu)化”平臺;第二、解決好FLUENT環(huán)境下水面船舶數值模擬過程中的不穩(wěn)定問題,實現(xiàn)“自動優(yōu)化”所要求的無人工干預目標。

        圖11 潛艇模型優(yōu)化前后阻力計算和試驗結果對比Fig.11 Computational and experimental resistance results of original and optimized hulls

        [1]楊佑宗,楊 奕,陳文煒.船舶線型設計與研究[J].上海造船,2001(2):18-23.

        [2]Yang C,Lohner R,Soto O.Optimization of a wave cancellation multihull ship using CFD tools[J].Journal of Hydrodynamics,2002,14(1):1-7.

        [3]Janson C E,Larsson L.A method for the optimization of ship hulls from a resistance point of view[R].Doktorsavhandlingar vid Chalmers Tekniska Hogskola,1997:35.

        [4]Larsson L,Baba E.Advances in Fluid Mechanics[M].Southampton,England:Computational Mechanics Publisher,1996(5):1-75.

        [5]Scott Percival,Dane Hendrix,Francis Noblesse.Hydrodynamic optimization of ship hull forms[J].Applied Ocean Research,2001,23(6):337-355.

        [6]Tahara Y,Stern F,Himeno Y.Computational fluid dynamics-based optimization of a surface combatant[J].Journal of Ship Research,2004,48(4):273-287.

        [7]Nancy C Groves,Thomas T Huang,Ming S Chang.Geometric characteristics of DARPA Suboff models[R].DTRC/SHD-1298-01,1989.

        [8]Han Lieh Liu,Thomas T Huang.Summary of DARPA Suboff experimental program data[R].Naval Surface Warfare Center,Carderock Division,1998:1-24.

        [9]傅慧萍,繆國平,高霄鵬.潛體繞流及遠場聲特性分析[J].海洋工程,2005,23(3):60-64.

        [10]劉文卿.實驗設計[M].北京:清華大學出版社,2005.

        [11]趙選民.實驗設計方法[M].北京:科學出版社,2006.

        [12]Noesis Solutions.Optimus Theoretical Background[M].Belgium:NOESIS SOLUTIONS,2006.

        [13]陽明盛,羅長童.最優(yōu)化原理、方法及求解軟件[M].北京:科學出版社,2006.

        Hull lines automatic optimization based on design space exploration

        LIANG Jun1,XU Jin-song1,XIE Jie1,YANG Xiao-yu2
        (1 School of Naval Architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,Shanghai 200030,China;2 Shanghai Karmon Technology Co.,Ltd.,Shanghai 200235,China)

        Ship hull optimization is the most important procedure to improve the resistance performance and operational economy.In this paper,an ‘Automatic Optimization’ platform for the ship hull optimization is presented by integrating the hull form modification,CFD computation,design space exploration,and design optimization together.This platform is applied to the bow optimization of the submarine DTRC Model 5470 from USA Naval Surface Warfare Center.The optimization results are validated through the model towing tests,and show the feasibility and effectiveness of the presented ‘Automatic Optimization’ platform.

        automatic optimization;ship hull form;resistance performance

        U661.31

        A

        1007-7294(2010)07-0741-08

        2008-12-29 修改日期:2009-12-22

        海洋工程國家重點實驗室自由研究課題資助(GKZD010012)

        梁 軍(1983-),男,上海交通大學船舶海洋與建筑工程學院碩士研究生。

        猜你喜歡
        型線潛艇阻力
        十分鐘讀懂潛艇史(下)
        鼻阻力測定在兒童OSA診療中的臨床作用
        潛艇哥別撞我
        十分鐘讀懂潛艇史(上)
        潛艇躍進之黃金時代
        零阻力
        英語文摘(2020年10期)2020-11-26 08:12:12
        高次曲線組合型線渦旋盤性能研究*
        機械制造(2020年8期)2020-09-30 06:32:24
        型線絞合導體ZC-YJLHV22-103×630鋁合金電纜的設計和生產
        電線電纜(2018年3期)2018-06-29 07:41:00
        別讓摩擦成為學習的阻力
        變截面復雜渦旋型線的加工幾何與力學仿真
        精品人妻一区二区三区蜜臀在线| 色狠狠av老熟女| 国产成人av一区二区三区无码| 免费 无码 国产精品| 亚洲日本中文字幕乱码| 日韩 无码 偷拍 中文字幕| 少妇高潮尖叫黑人激情在线| 亚洲两性视频一三区| 国产精品亚洲最新地址| 草草影院ccyy国产日本欧美| 成人免费看片又大又黄| 亚洲国产精品久久久久久网站| 中文字幕丰满人妻有码专区| 国产亚洲成人av一区| 男男受被攻做哭娇喘声视频| 在线视频一区二区日韩国产| 亚洲综合伊人久久综合| 国产精品无码一区二区三区在| 亚洲一区二区三区av资源| 国产成人无码综合亚洲日韩| 久久99精品国产99久久| 亚洲色偷偷偷综合网另类小说| av成人一区二区三区| 一本本月无码-| 又粗又大又黄又爽的免费视频| 精品免费久久久久国产一区| 成av人片一区二区久久| 四虎影视成人永久免费观看视频| 精品国产高清a毛片无毒不卡| 中文字幕日本一区二区在线观看| 国产精品网站91九色| 日韩毛片免费无码无毒视频观看| 亚洲人妻无缓冲av不卡| 最新国产精品国产三级国产av| 色欲欲www成人网站| 麻豆国产人妻欲求不满| 99精品国产av一区二区| 国产美女主播视频一二三区 | 欧美性狂猛xxxxx深喉| 91精品在线免费| 开心久久综合婷婷九月|