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

旋翼翼型氣動設計與評估軟件HRADesign

2021-09-17 08:15:34孫俊峰盧風順
空氣動力學學報 2021年4期
關鍵詞:優化方法模型

孫俊峰,盧風順,黃 勇,江 雄,牟 斌,許 勇

(中國空氣動力研究與發展中心,綿陽 621000)

0 引 言

直升機由于其獨特的飛行能力受到各航空大國的普遍重視。旋翼系統是直升機的關鍵部件,是直升機主要的升力面、推力面和操縱面。旋翼槳葉是由翼型構成的,翼型和槳葉氣動外形對旋翼性能有重要影響,高性能旋翼翼型可以提高旋翼懸停效率3%-5%,對直升機的前飛速度、等效升阻比和機動飛行能力、噪聲水平等都有很大影響。高性能旋翼翼型設計技術也一直是直升機設計的核心技術,是衡量直升機技術水平的重要標志。

在不同的飛行狀態和不同的槳葉半徑位置上,旋翼槳葉翼型的運行環境是迥然不同的[1]。翼型的設計需要在較寬的馬赫數范圍內,有較高的靜態和動態最大升力系數,以適應機動過載狀態;在較高的馬赫數及小迎角時,有較大的阻力發散馬赫數,以推遲前行槳葉激波失速;在中等馬赫數及中等迎角時,有較高的升阻比,以提高旋翼的懸停效率;在較低的馬赫數及大迎角時,有較好的失速特性,以延緩后行槳葉的氣流分離;在整個飛行包線內,有較小的俯仰力矩系數,以降低槳葉的操縱載荷,旋翼翼型設計需要綜合考慮不同的飛行環境,具有明顯的多點、多目標、強約束的特點。

旋翼翼型的設計方法主要分為反設計方法和數值優化設計方法兩大類。反設計方法主要是根據給定的速度分布或壓力分布設計翼型,缺點是理想的翼型速度分布、壓力分布特性難以預先給定。隨著CFD技術和計算機技術的不斷進步,數值優化設計類方法逐漸成為翼型設計的主流方法,該類方法可以直接實現多目標氣動性能的優化,降低了對人經驗的依賴,符合旋翼翼型設計多點/多目標的應用環境。數值類方法又包含控制論方法[2-3]和進化算法[4-5]兩大類,分別適用于不同的優化環境。

從20世紀70年代開始,直升機技術發達國家都開始致力于新型旋翼高性能翼型的研究,發展了一系列旋翼專用翼型,如美國Boeing-Vertol公司的VR翼型族、Sikorsky公司的SC翼型族,法國宇航院(ONERA)的OA翼型族,俄羅斯的TsAGI翼型族等[6-7]。這些翼型的成功應用,對改善直升機性能起到了重要作用。目前這些翼型族仍在繼續發展中,但鮮有相關文獻報道。國內也從20世紀90年代開始了對自主旋翼翼型的相關設計研究工作。尚克明等[8]基于Euler方程采用反設計的方法進行了旋翼翼型的設計研究。劉剛等[9]基于進化算法和Kriging模型開展了旋翼翼型多目標優化設計研究,優化結果經過風洞試驗驗證,滿足了設計要求。楊旭東等[10]基于梯度信息改進的響應面方法建立了旋翼翼型多點多約束氣動優化策略。韓忠華等[11]發展了基于Kriging模型與遺傳算法的旋翼翼型多目標多約束氣動優化設計方法,應用于OA209的設計取得了較好的效果。招啟軍等[12]開展了考慮旋翼翼型定常-非定常影響的綜合優化設計研究,優化結果在設計狀態下明顯改善了翼型的動態失速特性。

近年來直升機工業得到了長足的發展,先進直升機對高性能旋翼翼型的需求愈發迫切。目前國內旋翼翼型的研究正處于探索性研究向旋翼翼型系列研制的起步發展階段,自主翼型系列的研制還缺乏通用、高效、魯棒的設計工具,還沒有自主翼型應用于型號研究的先例。高性能旋翼翼型的設計,需要有先進的設計和優化框架作為支撐,以框架體系為基礎消化吸納各類先進技術,促進優化系統的不斷迭代更新完善,才能更有效地促進旋翼翼型工程實際設計能力的不斷提升。

本文通過綜合采用基于進化算法的多目標優化設計方法、基于主成分分析的多目標降維技術、翼型參數化技術以及高精度CFD性能分析工具,自主開發了旋翼翼型氣動設計與評估軟件系統——HRADesign。本文主要介紹了系統的架構設計和工作流程,詳細介紹了系統的多目標優化方法模塊、多目標降維技術、幾何管理模塊以及CFD性能分析模塊等功能模塊。利用該系統對ADODG標準翼型優化算例、旋翼翼型的常規多目標優化設計以及考慮多目標降維的旋翼翼型優化設計進行了研究,驗證了優化設計軟件系統的有效性。最后部分對本文進行了總結。

1 軟件架構設計

HRADesign作為一個綜合的旋翼翼型設計和評估的軟件系統,系統開發基于Eclipse集成環境,選用Python和C++語言作為開發工具,XML文件用于交換信息,采用wxWidgets軟件包提供可視化支持。用戶可以在PC機上完成優化問題的描述和參數的輸入后,通過集群系統完成優化設計的流程。系統根據高內聚、低耦合的原則分成層次結構,包含用戶界面層、功能模塊管理層、數據傳輸層、以及基礎服務層。圖1給出了軟件的體系結構示意圖,各層的功能定位如下:

圖1 軟件體系結構Fig.1 A diagram of the software architecture

1)用戶界面層。界面層給用戶提供了友好的人機交互接口來操作整個軟件,借助該層功能,用戶可實施優化設計問題建立、優化進度監控、優化結果查驗等操作。

2)功能模塊管理層。優化設計層集成了優化方法、代理模型、翼型參數化、優化目標處理、PCA降維技術、專家系統等模塊,借助CFD或試驗數據庫的數據支撐,完成各種氣動外形的優化設計功能。

3)數據傳輸層。該層主要負責為自適應優化設計層生產、收集和存儲數據,提供目錄管理、文件傳輸等基礎功能。

4)基礎服務層。該層提供整個軟件運行所需要的基礎服務,包括軟件服務(CFD解算器、網格程序、數據庫系統)和硬件服務(本地計算資源和大規模集群資源)。

圖2給出了HRADesign系統進行優化設計流程的示意圖。初始翼型進入系統,進行參數化并提取設計變量;選擇試驗設計方法,對每個樣本點進行網格生成和氣動特性分析;用氣動分析得到的數據來構建代理模型,利用該模型分析優化問題的優化目標和約束;最后進入優化迭代流程,直至獲得最終的優化翼型外形。

圖2 系統優化設計流程示意圖Fig.2 A diagram of the optimization process

2 系統功能模塊

HRADesign系統由多目標優化方法模塊、翼型參數化管理模塊、代理模型模塊、CFD性能分析模塊等主要的功能模塊構成,此外還包含問題定義,顯示監控等輔助模塊。系統提供了基于Windows的圖形操作界面,用戶在前臺通過操作界面完成優化問題的定義以及各項參數的設置,系統底層封裝了與后臺集群系統的連接以及信息的交換,可以實現本地計算與集群系統計算的無縫切換,用戶可以通過輸入/輸出系統監視優化進程,方便處理各類操作問題。圖3給出了系統的主要應用界面示意圖。

圖3 軟件應用界面示意圖Fig.3 A diagram of the graphical user interface

2.1 多目標優化方法

多目標優化方法模塊主要完成各種全局和局部數值優化方法的封裝,包含有進化算法和伴隨方法等,用戶通過應用界面可以選擇不同的優化方法完成翼型的魯棒設計、PCA降維分析等多目標設計功能。

HRADesign系統是以優化算法為核心,驅動整個優化流程的發展來進行翼型的優化設計。系統采用的多目標優化方法大多是以進化算法為基礎,結合多目標Pareto解的概念以及約束處理機制發展起來的,進化類算法屬于全局類優化方法,優化過程不依賴目標函數與設計變量的梯度信息,適合處理旋翼翼型復雜流動中的各類非線性問題。進化算法通過模擬生物種群的進化過程,利用選擇、交叉、變異等進化算子來找到多目標問題的Pareto解,解決了傳統優化設計中多目標加權的權重系數難以給定的難題,得到廣泛應用。圖4給出了利用該方法進行優化設計的流程圖。

圖4 進化算法流程圖Fig.4 The flow chart of the optimization process based on the evolutionary algorithm

2.2 主成分分析

主成分分析PCA的方法[13]是Deb在2005年提出的用于多目標降維的算法。算法首先對數據集進行標準化處理,求解目標函數的相關矩陣,得到相關矩陣的特征值及其特征向量,將特征值按從大到小的順序排序并求出其貢獻率,當貢獻率的累積大于初始給定的參數閾值時停止累積。閾值的選取在很大程度上影響實驗結果,如果閾值取太大,可能所有目標都會被選取,如果太小又容易丟失非冗余目標的信息。文獻一般建議閾值取值為95%。算法簡介如下[14]:

1) 設置閾值,非冗余目標集合Ⅰ為空。

2) 計算數據集。隨機初始化種群,對種群中的個體計算目標集中的所有目標,得到個體目標的數據集合P。

3) 對數據集P進行PCA處理,選取非冗余目標Ⅰ:

(a) 對目標集進行歸一化處理,計算相關矩陣,對相關矩陣計算特征值及對應的特征向量,依據特征向量來選取目標。

(b) 將特征值按照從大到小的順序依次排序,計算每個特征值占總特征值的比率,則每個特征值對應的特征向量依次稱為第一、第二、···、第n主成分。

(c) 依次分析每個主成分。如果主成分里的元素有正有負,選取最大和最小元素對應的目標加入非冗余目標集Ⅰ;如果所有元素均為正,選取最大元素值對應的目標加入非冗余目標集Ⅰ;如果所有元素均為負,將所有目標均加入到非冗余目標集Ⅰ。

(d) 如果主成分對應的特征值小于0.1,選取最大元素對應的目標加入到非冗余目標集Ⅰ。

(e) 考察特征值比率,如果大于閾值,停止分析過程,輸出非冗余目標解集Ⅰ;否則轉到(c)繼續非冗余目標的選取。

在優化設計中通過PCA主成分分析得到各目標之間的關聯關系后,提取決定問題本質的主要目標,將冗余目標剔除、或者轉化為約束條件,將高維多目標優化轉化為低維優化問題??梢越鉀Q多目標優化收斂慢甚至不收斂的問題,提高優化結果的可靠性。圖5給出了基于PCA分析的優化設計流程圖。

圖5 基于PCA分析的多目標優化設計流程圖Fig.5 The flow chart of PCA

DTLZ測試函數[15]是由Deb提出的一組測試多目標優化算法性能的測試函數,共有9組函數,DTLZ測試函數有已知的Pareto最優解,其中DTLZ5(I,M)(M表示目標個數,I表示非冗余目標個數),用于測試算法處理包含冗余目標的能力。選取DTLZ5(2,10)函數用于測試,該函數有10個目標,其中2個為非冗余目標,函數公式見式(1):

對該函數進行優化測試,經過PCA分析,得到目標9和目標10兩個非冗余目標。首先對10個目標進行優化,得到優化結果對應的兩個非冗余目標的Pareto前沿如圖2(a)所示,去除冗余目標后再進行優化,得到優化結果如圖2(b)所示。DTLZ5(2,10)函數的Pareto前沿收斂到圓弧曲線[16],通過圖6可以看出,算法在處理冗余目標后,達到了最優Pareto前沿。

圖6 DTLZ5(2,10)測試函數收斂比較Fig.6 Convergence of the test function DTLZ5(2,10)

2.3 翼型參數化

翼型參數化管理模塊主要完成翼型的輸入/輸出,約束評估、翼型曲線參數化和網格自動重構、設計變量的選取和確定設計空間范圍等功能。目前常用的參數化方法有解析函數線性疊加法,NURBS曲線[17]、CST[18]方法、FFD[19]方法等。本系統主要采用CST技術實現翼型的參數化表示。

CST方法是波音公司B.M.Kulfan等提出的一種通用幾何參數化表示方法??梢杂媒y一的解釋函數表示鈍前緣/尖后緣類翼型和雙鈍頭翼型等新型翼型。用該方法描述翼型,容易控制前緣半徑、彎度/厚度分布、后緣角以及后緣厚度等關鍵參數,而且設計參數數目容易控制,也具有局部修改控制的能力。

翼型上下表面均用以下公式描述:

其中類型函數:

形狀函數:

翼型表面形狀改變后采用基于雙曲方程[20]的網格生成方法完成計算網格重構,可以直接輸出給CFD解算器使用。

2.4 Kriging代理模型

優化設計過程中,采用高精度CFD分析工具進行翼型氣動性能計算、優化目標函數和約束函數的評估,會帶來計算資源成本過高、計算周期過長的問題。為了提高設計的效率,可以采用代理模型的方法構造目標函數和約束函數的近似函數,優化算法作用于近似目標函數和近似約束函數以尋找優化問題的最優解,通過近似函數的不斷改進和優化算法的不斷迭代,直至最終滿足收斂條件。

HRADesign系統中,代理模型技術主要包括試驗設計方法和代理模型方法兩部分。試驗設計方法決定了樣本點的個數和樣本點的空間分布情況,系統采用拉丁超立方采樣[21]和均勻采樣[22]的方法選取樣本點,以保證樣本點在設計空間的均勻分布。Kriging模型[23-24]則用來作為目標函數和約束函數的近似模型。

Kriging代理模型起源于地理空間統計學,是一種估計方差最小的無偏估計插值模型,具有全局近似和局部隨機誤差估計相結合的特點。通過Kriging模型可以得到未知點的函數值和不確定性,因而在優化過程中,Kriging模型需要根據優化進程自適應地更新,增加樣本點在非線性區域的分布,提高模型的預測精度,系統采用EI方法或者最小值加點準則等方法,提高模型的精度和自適應能力。圖7給出了基于代理模型的優化流程示意圖。

圖7 基于代理模型的優化流程示意圖Fig.7 The flow chart of a surrogate-based multi-objective optimization

第一步:試驗設計。采用均勻設計等試驗設計方法在設計空間中選取樣本點。

第二步:樣本點性能評估。對樣本點分別生成計算網格,采用CFD工具進行氣動性能評估。

第三步:構建Kriging模型。利用樣本點的性能計算結果構建初始Kriging代理模型。

第四步:利用構建的Kriging代理模型,結合基于進化算法的多目標優化方法進行設計優化流程。

第五步:Kriging模型的重構。優化過程中評估種群個體的EI值,根據種群個體的EI值選定附加的樣本點,對樣本點進行計算網格的生成和CFD性能評估,在N=N+ 1個樣本點的基礎上進行Kriging模型的重構。

第六步:返回第四步繼續迭代循環。

第七步:直到系統收斂或達到設定的優化步數,返回優化結果。整個優化流程結束。

2.5 氣動性能分析

旋翼翼型CFD性能計算的精度和效率是進行氣動設計優化的關鍵,對優化設計的結果有直接影響。

HRADesign系統中采用的CFD分析工具是自主開發的二維RANS解算器MBNS2D。MBNS2D采用格心型有限體積方法求解雷諾平均Navier-Stokes方程,空間離散采用Roe格式,湍流模型包含SA一方程模型和SST兩方程模型,通過多重網格技術進行流場加速收斂,提高計算效率。

在旋翼翼型的計算中,轉捩對翼型前緣流動影響很大,如果不考慮轉捩的影響,阻力及零升阻力的計算結果會與試驗差別很大,為了增強方法的適應性,HRADesign系統中采用了 γ -Reθ湍流轉捩模型[25]來提高翼型阻力系數的計算精度。

3 優化算例及分析

3.1 算例1

為了提高優化設計軟件的可信度,AIAA氣動優化設計討論組(ADODG)給出了一套優化設計的標準算例,用于優化設計軟件的驗證和確認。這些算例包含了翼型和機翼在氣動和幾何約束條件下阻力最小化優化設計問題。本文以其中跨聲速條件下RAE2822翼型阻力最小化算例來考核優化設計系統的性能。

以RAE2822翼型作為初始翼型開始優化,翼型的跨聲速設計條件為:Ma∞= 0.734,Re= 6.5×1 06,CL= 0.824,選取阻力系數CD最小化作為優化目標,幾何約束條件為保持優化后翼型的面積不減少。優化問題數學描述為:

圖8和圖9分別給出了優化前后翼型的外形和表面壓力分布的比較,圖10給出了優化前后翼型壓力云圖的比較??梢钥闯?,優化翼型前緣吸力峰增強,基本消除了初始翼型中段的強激波,等值線變得平順光滑。表1給出了翼型在優化前后的氣動特性比較,在保持升力和力矩的氣動約束以及面積不減的幾何約束條件下,優化翼型的阻力系數降低了約87 counts,優化翼型的阻力特性得到明顯改善。文獻[26]對該算例給出了進一步的分析。

表1 優化前后翼型性能比較Table 1 Performance of the original and optimized RAE2822 airfoils

圖8 優化前后翼型的形狀比較Fig.8 Configurations of the initial and optimized RAE2822 airfoils

圖9 優化前后翼型表面壓力分布比較Fig.9 The comparison of airfoil pressure distributions between the initial and optimized airfoils

圖10 優化前后翼型壓力云圖比較Fig.10 The comparison of pressure contours between the initial and optimized airfoils

3.2 算例2

直升機旋翼翼型的設計,需要考慮多種飛行狀態下的性能,一般要求:較高的最大升力系數;較高的阻力發散馬赫數;在較大的馬赫數范圍內有較高的升阻比以及較小的俯仰力矩,屬于典型的多目標設計優化問題。以某翼型為基本翼型,在保持厚度不減的條件下開展多點設計優化,使得前飛、機動、懸停條件下滿足設計指標的要求。優化問題描述為:

其中Maddm為阻力發散馬赫數,K為升阻比,下標0表示基本翼型的性能指標。

圖11給出了優化翼型和基本翼型外形的比較,優化外形頭部半徑增大,上表面厚度增加。

圖11 優化前后翼型的形狀比較Fig.11 Configurations of the initial and optimized airfoils

圖12給出了機動狀態Ma= 0.4下升力特性的比較,可以看出優化翼型最大升力系數有明顯提高。圖13給出了懸停狀態下極曲線的對比,可以看出優化翼型的阻力系數有所降低。表2給出了優化翼型和基本翼型的性能結果比較,從結果可以看出,優化翼型的機動和懸停性能有明顯提升,前飛性能在保持力矩性能的前提下阻力發散馬赫數略有降低。該厚度翼型主要位于槳葉的中段,在滿足厚度和力矩約束的條件下,需要兼顧低速和高速性能,要提高前飛阻力發散馬赫數,就要損失低速升力性能,所以多目標優化設計中往往需要在性能之間進行折衷平衡。

表2 優化前后翼型性能比較Table 2 Performance of the initial and optimized airfoils

圖12 翼型升力特性比較( M a=0.4, R e=2.8×106)Fig.12 The comparison of lift coefficients between the initial and optimized airfoils ( M a=0.4,Re=2.8×106)

圖13 優化翼型與初始翼型極曲線特性比較(Ma=0.6, Re=4.2×106)Fig.13 The comparison of polars characteristics between the initial and optimized airfoils ( M a=0.6, R e=4.2×106)

3.3 算例3

以某12%厚度翼型作為參考的基本翼型,綜合利用本文發展的各項技術,開展多目標/多點設計優化,利用PCA方法實現多目標降維,驗證了多目標降維技術的有效性。問題可描述為:

上述優化問題共有12個目標,其中前飛狀態有2個目標,Mdd和Cm;機動狀態共有8個,分別對應4個馬赫數下的CLmax和Kmax;懸停狀態有2個目標,對應兩個狀態下的升阻比。

利用拉丁超立方采樣在設計空間中隨機抽樣,本文選取了480個樣本點,利用CFD方法對樣本點分別進行12個目標的氣動性能計算,得到480個樣本點的目標性能集,對性能集進行PCA分析。根據PCA的目標降維選取方法,對原始目標集進行降維后得到新的目標集為:

經過PCA分析可以看出,原始優化問題含有7個冗余目標,非冗余優化目標剩下5個,簡化了原始設計問題,針對經PCA降維分析后的問題開展研究。

圖14給出了優化翼型和基本翼型的外形比較,表3給出了兩者性能指標的比較。從結果可以看出,優化翼型略微變薄,最大厚度位置前移,優化翼型的前飛性能和懸停性能優于基準翼型,綜合性能較基本翼型有所提高,可以有效提高旋翼的懸停和機動性能。

圖14 初始翼型和優化翼型外形比較Fig.14 Configurations of the initial and optimized airfoils

表3 基準翼型和優化翼型性能比較Table 3 The performance of the initial and optimized airfoils

4 結 論

HRADesign軟件系統主要提供通用、魯棒、高效的直升機旋翼翼型工業設計能力。通過軟件工程的設計方法,構建了通用的旋翼翼型氣動設計和評估軟件系統。

1)系統采用軟件工程方法實現了架構設計和功能模塊的集成,通過分層設計實現了功能模塊的高內聚和低耦合,實現了功能模塊的靈活擴充,提高了系統可擴展性和用戶解決問題的靈活性,滿足了設計目標要求。

2)系統集成了多目標進化算法、CFD性能分析工具、翼型參數化工具、網格自動重構、Kriging代理模型以及多目標降維等功能模塊,可以滿足旋翼翼型多目標優化設計的功能需求。

3)考核算例及應用驗證了HRADesign集成系統功能模塊的有效性和通用、魯棒、高效的設計能力。PCA方法的應用展現了其在高維氣動多目標優化問題中的應用潛力,對于分析問題的主要特征,降低設計的復雜程度有重要指導意義。

4)HRADesign系統平臺目前主要集中于旋翼翼型靜態氣動特性多目標優化設計,下一步要發展完善靜態/動態特性多目標優化設計方法,進一步提高旋翼翼型的綜合設計能力。

猜你喜歡
優化方法模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
主站蜘蛛池模板: 六月婷婷精品视频在线观看| 亚洲愉拍一区二区精品| 亚洲IV视频免费在线光看| 伊人久久久久久久久久| 尤物特级无码毛片免费| 国产亚洲高清视频| 欧美精品亚洲精品日韩专区va| 中美日韩在线网免费毛片视频| 欧美午夜在线视频| 亚洲第一极品精品无码| 亚洲无码免费黄色网址| 久久国产精品嫖妓| 国产尤物jk自慰制服喷水| 久久精品电影| 国产在线观看人成激情视频| 老色鬼久久亚洲AV综合| 4虎影视国产在线观看精品| 国产精品自在自线免费观看| 亚洲天堂啪啪| 在线无码av一区二区三区| 狠狠色香婷婷久久亚洲精品| 99久久国产精品无码| 麻豆国产在线观看一区二区| 精品天海翼一区二区| 国产精品一线天| 欧美成人精品在线| 丝袜无码一区二区三区| 国产区在线观看视频| 久久鸭综合久久国产| 国产肉感大码AV无码| 亚洲欧美国产视频| 国产精品网址在线观看你懂的| 久草网视频在线| 亚洲精品不卡午夜精品| 青青久视频| 久久青青草原亚洲av无码| 欧美激情第一欧美在线| 欧美精品黑人粗大| 素人激情视频福利| 国产成人无码综合亚洲日韩不卡| 99视频全部免费| 日本精品一在线观看视频| 亚洲免费毛片| 欧美中出一区二区| 亚洲资源站av无码网址| 久操线在视频在线观看| 国产最新无码专区在线| 国产偷倩视频| 亚洲一区第一页| 午夜色综合| 成人综合网址| 日韩AV无码一区| 国产主播在线观看| 91在线国内在线播放老师| 免费在线播放毛片| 亚洲av日韩综合一区尤物| 久久综合九九亚洲一区| 久久伊人色| 一区二区午夜| 亚洲无码视频一区二区三区 | 青青热久免费精品视频6| 99久久精品国产麻豆婷婷| 青草国产在线视频| 99久久亚洲综合精品TS| 九九热免费在线视频| 日韩免费毛片| 在线观看视频99| 久久黄色视频影| 久久综合成人| 久久中文电影| 超清人妻系列无码专区| 久久福利网| 日韩在线观看网站| 亚洲av色吊丝无码| 91视频首页| 久久婷婷六月| 国产成a人片在线播放| 欧美精品成人| 一本大道香蕉久中文在线播放| 亚洲欧洲天堂色AV| 91美女视频在线观看| 在线观看亚洲国产|