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

面向乘波體氣動外形設(shè)計的代理模型特性分析及優(yōu)化

2015-04-22 05:50:58季廷煒王崯矚張繼發(fā)
固體火箭技術(shù) 2015年4期
關(guān)鍵詞:優(yōu)化模型設(shè)計

季廷煒,王崯矚,張繼發(fā),張 帥

(1.浙江大學 航空航天學院,杭州 310027;2. 浙江大學 工程與科學計算中心,杭州 310027)

?

面向乘波體氣動外形設(shè)計的代理模型特性分析及優(yōu)化

季廷煒1,2,王崯矚1,2,張繼發(fā)1,2,張 帥1,2

(1.浙江大學 航空航天學院,杭州 310027;2. 浙江大學 工程與科學計算中心,杭州 310027)

代理模型在氣動布局優(yōu)化設(shè)計得到了廣泛應用,研究探討分析各代理模型的特性分析是重要的基礎(chǔ)問題。文中基于3種Kriging和二階多項式響應面模型等4種代理模型,以升阻比和容積率為樣本,研究分析了其在乘波體優(yōu)化設(shè)計中的的擬合精度和結(jié)構(gòu)特性。研究表明,在高超乘波體優(yōu)化設(shè)計中,基于自然指數(shù)型評估函數(shù)的Kriging模型具有相對突出的擬合精度。最后,在乘波體升阻力系數(shù)和容積率代理模型的多目標優(yōu)化中,對比了基于自然指數(shù)型評估函數(shù)的Kriging模型和常用的基于高斯函數(shù)的Kriging模型的優(yōu)化結(jié)果,發(fā)現(xiàn)基于自然指數(shù)型評估函數(shù)的Kriging模型對于峰值的擬合效果更好。

代理模型;氣動布局;乘波體構(gòu)型;飛行器概念設(shè)計;多目標優(yōu)化

0 引言

高超聲速飛行器乘波體布局優(yōu)化設(shè)計過程中,存在高維性、多峰性、全局性和仿真分析計算量大等問題。為了高效準確地進行高超聲速氣動布局優(yōu)化,將代理模型方法引入高超聲速飛行器布局優(yōu)化設(shè)計,減少計算機高強度的仿真計算次數(shù),降低設(shè)計空間搜索的計算量,消除數(shù)值噪聲,提高全局尋優(yōu)效率[1-5]。

當前,有關(guān)氣動布局代理模型的分析和比較大部分針對機翼外形、飛行器和火箭等工程優(yōu)化上的應用,其中所提出的代理模型的比較研究大多集中于二階多項式響應面和Kriging代理模型、Kriging模型同多項式響應面模型和徑向基RBF插值模型之間的比較研究[6-14]。在高超飛行器工程應用方面,Tsuchiya[15]、車競[16]和陳萬春等[17]國內(nèi)外學者[18-27]都用代理模型構(gòu)建了相應的學科的響應面模型進行優(yōu)化設(shè)計,做出了突出的成果。其對應的研究主要針對響應面模型在高超聲速飛行器優(yōu)化中的構(gòu)建,亟需深入研究適用于高超聲速飛行器優(yōu)化設(shè)計的代理模型特性。

目前,文獻中高超聲速飛行器優(yōu)化設(shè)計中涉及的Kriging模型大多采用高斯模型(Gaussian)作為相關(guān)性評價函數(shù),對于其他相關(guān)性評價函數(shù)的模型研究的開展較少。本文以乘波體升阻力和容積為樣本,分別編寫了二階多項式響應面模型(QRS)、基于高斯相關(guān)評估函數(shù)的Kriging模型(GKG)、基于自然指數(shù)相關(guān)評估函數(shù)的Kriging模型(EKG)、基于一般指數(shù)相關(guān)評估函數(shù)的Kriging模型(XKG)等4種代理模型,探索適于乘波體氣動布局概念設(shè)計的代理模型,分析對比各設(shè)計參數(shù)對設(shè)計響應的影響情況。最后,采用遺傳多目標優(yōu)化算法,以升阻比和容積率作為優(yōu)化目標,應用代理模型技術(shù)解決高超聲速乘波體氣動布局優(yōu)化設(shè)計。

1 基本理論及方法

1.1 代理模型方法

目前,工程上常用的代理模型有二階多項式響應面模型和Kriging模型。

(1)二階多項式響應面模型

多項式響應面模型利用樣本數(shù)據(jù)擬合得出設(shè)計變量與性能響應間的多項式函數(shù)關(guān)系。多項式響應面計算量小,基本上能反映設(shè)計空間的非線性,一般采用二階多項式響應面模型,即

(1)

式中y為響應;x為設(shè)計變量;β為多項式參數(shù);ε為擬合誤差;n為設(shè)計變量的維數(shù)。

(2)Kriging模型

Kriging模型是來源于地質(zhì)學,結(jié)合了工程和統(tǒng)計的一種函數(shù)近似方法,具有相當?shù)撵`活性,但計算量較大。模型如下:

y(x)=f(x)+Z(x)

(2)

式中y(x)為待定響應;f(x)為多項式函數(shù),是對設(shè)計空間的全局近似;Z(x)是均值為0、方差為σ2、協(xié)方差非0的隨機過程,表示對全局近似的背離。

(3)

式中r(x)為輸入x與各樣本點{x,x2,…,xn}的相關(guān)系數(shù);R為樣本點相關(guān)矩陣。

r和R的基本元素為相關(guān)評估函數(shù)。最常用的是高斯相關(guān)評估函數(shù)。此外,還有自然指數(shù)相關(guān)評估函數(shù)和一般指數(shù)相關(guān)評估函數(shù)。

1.2 氣動性能評估方法

由于雙橢球模型在高超聲速流動中呈現(xiàn)復雜的流場特性,存在三維弓形頭激波、二次激波以及三維分離流動、分離激波,是一種具有代表性的高超聲速算例,適于數(shù)值模擬程序的可靠性驗證。

本文利用自研的CFD程序,對文獻[28]提供的風洞實驗雙橢球標模開展數(shù)值計算驗證。雙橢球模型幾何外形方程[28]如下:

(1)水平橢球

(2)垂直橢球

(3)上半柱形

(4)下半柱形

基座橢球在x、y、z方向上的半軸長分別為157.90、39.47、65.79 mm,兩橢球中心重合,坐標原點設(shè)在橢球中心點O。

利用該實驗工況馬赫4.94、來流溫度370 K作為計算來流條件,選用SSTκ-ε模型,對流項離散采用三階MUSCL格式。圖1(a)、(b)將CFD對稱面的上表面和下表面的壓力分布與文獻的實驗結(jié)果進行了對比。可見,壓力變化與實驗值接近,說明了自研程序的可靠性。

2 基于代理模型的優(yōu)化設(shè)計

2.1 系統(tǒng)框架

基于常用的氣動布局優(yōu)化設(shè)計代理模型構(gòu)建流程的基礎(chǔ)上,本文將基于CFD氣動特性分析,構(gòu)造基于代理模型的氣動布局優(yōu)化流程,如圖2所示。具體步驟如下:

(1)使用拉丁超立方抽樣(Latin Hypercube Sampling, LHS)的實驗設(shè)計策略產(chǎn)生一系列設(shè)計變量樣本點;

(2)選擇部分設(shè)計樣本點,進行CFD計算,獲得CFD計算結(jié)果;

(3)根據(jù)CFD計算結(jié)果,得到輸入同輸出響應集,構(gòu)建對應的代理模型,并對該近似模型進行可信度驗證后,執(zhí)行優(yōu)化設(shè)計;

(4)獲得氣動布局優(yōu)化結(jié)果后,利用CFD對優(yōu)化結(jié)果進行驗證計算,如結(jié)果不能滿足精度要求,將其結(jié)果用于樣本集合和代理模型的優(yōu)化修正;

(5)驗證最優(yōu)化點的精度,獲得合適的外形氣動布局結(jié)果。

(a)上表面壓力

(b)下表面壓力

2.2 實驗設(shè)計

乘波體的設(shè)計過程詳細生成過程見文獻[2, 4],設(shè)計馬赫數(shù)Ma、長寬比LWR、激波半錐角θ、底部曲線夾角ψ以及二次項系數(shù)B,確定乘波體外形。不失一般性,本文將馬赫數(shù)和長寬比設(shè)為固定值,分別為Ma=6和LWR=3,其他的初始設(shè)計參數(shù)為θ=3.948°,ψ=2.035°,B=0.771。為了獲得覆蓋較大的設(shè)計空間,將激波半錐角、底部曲線夾角、二次項系數(shù)將從區(qū)間值θ、ψ、B分別設(shè)為[0.05°,5°]、[0.05°,5°]、[0.05,0.5]。

為了探索樣本集數(shù)對代理模型建立的影響,利用LHS方法[19],分別選擇了樣本數(shù)31、64、91、126生成樣本點,構(gòu)建乘波外形,進行升阻力系數(shù)和容積評估。

CFD計算時,將飛行高度H=25 km的大氣環(huán)境作為計算來流的初始條件,選擇理想氣體模型,靜壓為2 549.2 Pa,靜溫為221.5 K,馬赫數(shù)為6。根據(jù)幾何外形左右對稱,采用一半計算域,基于結(jié)構(gòu)網(wǎng)格,在靠近乘波體表面對網(wǎng)格進行加密。法向壁面網(wǎng)格y+在1左右。圖3給出論文討論的某典型乘波體外形及計算網(wǎng)格,其他計算網(wǎng)格與其類似。

圖2 基于代理模型的優(yōu)化設(shè)計流程Fig.2 Flowchart of surrogate based on the design optimization

(a)乘波體外形 (b)計算網(wǎng)格

3 代理模型評估

描述代理模型精度的方法有很多,本文首先采用均方根誤差(Root Mean Square Error,RMSE):

(4)

此外,為了描述代理模型和氣動評估結(jié)果之間的相近程度,在本文利用皮爾氏相關(guān)系數(shù)r-correlation用于該項評估。

(5)

其中

(6)

(7)

(8)

3.1 升阻比系數(shù)代理模型分析

3.1.1 升阻比代理模型精度和相關(guān)性評估

圖4描述了樣本空間中隨著樣本數(shù)的增加,升阻比代理模型的RMSE和r-correlation值的變化情況。其中,水平軸表示取樣點數(shù)量,縱軸分別表示RMSE和r-correlation值。

由圖4(a)可知,當樣本數(shù)少時,XKG的RMSE值在這4種模型中最低,表示其擬合精度最好,而GKG效果最差;當樣本點數(shù)量達到126時,EKG是4種代理模型中的RMSE值最接近0,表示其對實際值的擬合精度最高。

由圖4(b)可知,與RMSE誤差值類似,當樣本點最少時,GKG模型的相關(guān)性評估是最差,而EKG模型的相關(guān)性最好;當樣本點為126時,所構(gòu)建代理模型r-correlation值最接近1的是EKG,表示了該模型和氣動評估結(jié)果之間的相近程度越好。

由表2可看出,GKG、XKG和EKG的相關(guān)性函數(shù)系數(shù)沒有隨樣本數(shù)的增加而明顯變化,說明樣本數(shù)的增加并沒有對代理模型的結(jié)構(gòu)產(chǎn)生很大的影響。

總之,針對升阻比系數(shù)的代理模型,根據(jù)RMSE值和r-correlation的分析,EKG相比與文獻中常用的GKG代理模型具有更好的擬合效果,在同樣的樣本集中能夠產(chǎn)生更好的擬合精度。

3.1.2 升阻比代理模型影響因素分析

表3列出了QRS隨樣本增加其各項系數(shù)的變化情況。可明顯看出,激波半錐角(β2)系數(shù)值最大,說明升阻比代理模型受激波半錐角(β2)影響最大,次之為底部曲線的二次項系數(shù)。

為了驗證其存在的影響程度,將EKG代理模型的預測值由圖5所示。可明顯發(fā)現(xiàn),隨著激波半錐角的變化,其多峰性明顯,表明乘波體升阻性能受激波半錐角的影響波動性強。

總之,由表3和圖5可共同發(fā)現(xiàn),乘波體升阻比代理模型的性能明顯受激波半錐角的影響,對于該類飛行器的設(shè)計有重要影響。

(a)RMSE

(b)r-correlation

表2 升阻比Kriging代理模型系數(shù)Table 2 Values of Kriging parameters for lift-to-drag ratio

表3 升阻比二次多項式響應面模型系數(shù)Table3 Values of quadratic polynomial regression parameters for lift-to-drag ratio

(a)31 樣本 (b)64樣本

(c)91樣本 (d)126樣本

3.2 容積率代理模型

3.2.1 容積率代理模型精度和相關(guān)性評估

圖6描述隨著樣本集合中樣本數(shù)的增加,容積代理模型的RMSE和r-correlation值的變化,同樣的水平軸表示樣本數(shù),縱軸表示RMSE和r-correlation值。同升阻比代理模型一樣,兩者的精度是比較接近的,即代理模型用于升阻比和容積兩者之間的精確程度類似,兩者之間的變化影響程度相當。

由圖6(a)可知,同樣可發(fā)現(xiàn),當樣本點數(shù)量少時,GKG性能最差,而EKG和XKG性能較好,RMSE誤差值下降很快,但存在波動比較大的問題。同樣的當樣本點數(shù)量達到126時, EKG是4種代理模型中的RMSE值最接近0,其對實際值的擬合精度最高。

如圖6(b)所示,當樣本點少時,QRS代理模型相比于其他3類具有明顯優(yōu)勢,當樣本數(shù)量增加時,EKG和XKG性能提高很快,但同樣存在波動較大的問題。當樣本點數(shù)量達到126時,EKG相比與其他代理模型具有更好的計算精度,與氣動評估值相關(guān)性最好。

由表4的系數(shù)變化情況,同樣可說明樣本數(shù)的增加并沒有對代理模型的結(jié)構(gòu)產(chǎn)生很大的影響。

(a)RMSE

(b)r-correlation

總之,針對容積率代理模型,同樣地從RMSE值和r-correlation來看,EKG相比與其他代理模型具有更好的擬合精度。

3.2.2 容積率代理模型影響因素分析

由表5列出了QRS隨樣本數(shù)增加其各項系數(shù)的變化情況。可明顯看出,激波半錐角(β2)系數(shù)值最大,說明容積率代理模型同樣受激波半錐角(β2)影響。

圖7描述的是GKG代理模型的預測值的變化情況。與升阻比代理模型相比,容積率代理模型受設(shè)計參數(shù)的影響程度上少一個數(shù)量級,但同樣可發(fā)現(xiàn),其激波半錐角峰值變化劇烈,表明乘波體容積同樣容易受激波半錐角的影響,其對升阻比影響程度更大。

表4 容積率Kriging代理模型系數(shù)Table4 Values of Kriging parameters for volumetric efficiency

表5 容積率二次多項式響應面模型系數(shù)Table5 Values of quadratic polynomial regression parameters for volumetric efficiency

(a)31樣本 (b)64樣本

(c)91樣本 (d)126樣本

4 優(yōu)化結(jié)果和討論

本文基于上述代理模型的研究結(jié)果,采用126個樣本點所構(gòu)建的升阻比、容積率的EKG和GKG代理模型分別作為優(yōu)化設(shè)計的目標函數(shù),選用多目標遺傳算法進行優(yōu)化,優(yōu)化結(jié)果如表6所示。

采用不同的代理模型作為優(yōu)化目標,可得到當126個樣本點生成的EKG模型作為優(yōu)化目標時,優(yōu)化升阻比和容積分別為3.122 6與0.848,相比于初始設(shè)計結(jié)果升阻比2.669和0.278,升阻比提高了16.99%,容積率提高了205.03%;當126個樣本點生成的GKG模型作為優(yōu)化目標時,優(yōu)化升阻比和容積分別為2.998 3與0.748, 相比于初始設(shè)計值,升阻比提高了12.33%,容積率提高了127.33%。

在同樣的構(gòu)建樣本下,采用一樣的優(yōu)化算法,初始參數(shù)也相同的情況下,EKG模型得出優(yōu)化結(jié)果比GKG的好。由此也可說明,在乘波體優(yōu)化設(shè)計中,EKG代理模型對峰值的擬合上,相比于文獻中常用的GKG模型精度上更高。

表6 優(yōu)化設(shè)計結(jié)果Table6 Optimization results

5 結(jié)論

(1)在乘波體優(yōu)化設(shè)計過程中,EKG相比于QRS、GKG和XKG的擬合精度較高,能夠產(chǎn)生了更佳的擬合效果。

(2)隨著樣本數(shù)的增加,4種代理模型的擬合精度也越好,但其對乘波體升阻力和容積率代理模型結(jié)構(gòu)影響較小。

(3)基于本文研究結(jié)論,采用EKG模型,實現(xiàn)基于代理模型的多目標氣動布局優(yōu)化方法,優(yōu)化結(jié)果與初始設(shè)計升阻比提高16.99%,容積率提高了205%。

[1] 羅世彬, 羅文彩, 王振國. 基于試驗設(shè)計和響應面近似的高超聲速巡航飛行器多學科設(shè)計優(yōu)化[J]. 導彈與航天運載計算, 2003, 266(6): 2-9.

[2] 王爍,李萍,陳萬春. 錐導乘波體氣動代理模型建模方法研究[J]. 飛行力學, 2012, 30(1): 43-47.

[3] 張珍銘, 丁運亮, 劉毅. 升力體外形設(shè)計的代理模型優(yōu)化方法[J]. 宇航學報, 2011, 32(7): 1435-1444.

[4] 耿永兵, 劉宏, 姚文秀,等. 錐形流乘波體優(yōu)化設(shè)計研究[J]. 航空學報, 2006, 27(1): 23-28.

[5] 許少華, 侯中喜, 葛愛學,等. 錐導乘波構(gòu)型設(shè)計、優(yōu)化與分析[J]. 推進技術(shù), 2008, 29(4): 448-453.

[6] 穆雪峰, 姚衛(wèi)星, 余雄慶,等. 多學科設(shè)計優(yōu)化中常用代理模型的研究[J]. 計算力學學報, 2005, 22(5): 608-612.

[7] Etaman L F P. Design and analysis of computer experiments: the method of Sacks et al.[R]. Engineering Mechanics report WFW 94.098. Eindhoven University of Technology, 1994.

[8] Giunta A A. Aircraft multidisciplinary design optimization using design of experiments theory and response surface modeling methods[D]. Virginia Polytechnic Institute and State University, 1997.

[9] Jin R, Chin W, Simpson T W. Comparative studies of metamodeling techniques under multiple modeling criteria[R]. AIAA 2000-4801.

[10] Simpson T W, Mauery T M, Korete J J, et al. Kriging models for global approximation in simulation-based multidisciplinary design optimization[J]. AIAA Journal, 2001, 39(12): 2233-2241.

[11] Simpson T W, Peplinski J D, Loch P N, et al. Metamodels for computer based engineering design: Survey and recommdation[J]. Engineering with Computers, 2001, 17(2): 129-150.

[12] Viana F A C, Simpson T W, Balabanov V, et al. Metamodeling in multidisciplinary design optimization: How far have we really come?[J]. AIAA Journal, 2014, 52(4): 670-690.

[13] Chung H S, Alonso J J. Comparison of approximation models with merit function for design optimization[R]. AIAA 2000-4754.

[14] Zhang K, Han Z, Li W, et al. Coupled aerodynamic/structural optimization of a subsonic-transport wing using surrogate model[J]. Journal of Aircraft, 2008, 45(6): 2167-2171.

[15] Tsuchiya T, Takenaka Y, Taguchi H. Multidisciplinary design optimization for hypersonic experimental vehicle[J]. AIAA Journal, 2007, 45(7): 1655-1662.

[16] 車競, 王文正, 何開鋒. 高超聲速飛行器并行子空間優(yōu)化設(shè)計[J]. 航空動力學報, 2010, 25(4): 912-917.

[17] 馬樹微, 杜斌, 陳萬春. 基于代理模型的高超聲速飛行器協(xié)同優(yōu)化設(shè)計 [J]. 北京航空航天大學學報, 2013, 39(8): 1042-1047.

[18] Jeong S, Murayama M, Yamamoto K. Efficient optimization design method using Kriging model[J]. Journal of Aircraft, 2005 42(2): 413-420.

[19] Lophaven S N, Nielsen H B, Sondergaard J. DACE: A MATLAB Kriging Toolbox[R]. Denmark University of Technology, 2002.

[20] Senant N E, Bloor M I G, Wilson M J. Aerodynamic design of a flying wing using response surface methodology[J]. Journal of Aircraft, 2000, 37(4): 562-569.

[21] Vavalle A, Qin N. Iterative response surface based optimization scheme for transonic airfoil design[J]. Journal of Aircraft, 2007, 44(2): 365-376.

[22] Yim J W, Lee B J, Kim C, et al. Multi-stage aerodynamic design of multi-body geometries via global and local optimization methods[R]. AIAA 2008-134.

[23] 鄧艷丹, 黃生洪, 楊基明,等. 一種X-51A相似飛行器模型的氣動特性初探[J]. 空氣動力學學報, 2013, 31(3): 376-380.

[24] Paiva R M, Carvalho A R D, Suleman A, et al. Comparison of Surrogate Models in a Multidisciplinary Optimization Framework for Wing Design[J]. AIAA Journal, 2010, 48(5): 995-1006.

[25] Koziel S, Leifsson L. Surrogate-based aerodynamic shape optimization by variable-resolution models[J]. AIAA Journal, 2013, 51(1): 94-106.

[26] Pourbagian M, Habashi W G. Surrogate-based optimization of electro-thermal wing anti-icing systems[J]. Journal of Aircraft, 2013, 50(5): 1555-1563.

[27] Ursache N M, Bressloft N W, Keane A J. Aircraft Roll Enhancement via Multi-Objective Optimization Using Surrogate Modeling[J]. AIAA Journal, 2011, 49(7):1525-1541.

[28] 李素循. 典型外形高超聲速流動特性[M]. 北京:國防工業(yè)出版社, 2008.

(編輯:呂耀輝)

Surrogate models characteristic analysis and optimization for waverider aerodynamic design

JI Ting-wei1,2,WANG Yin-zhu1,2,ZHANG Ji-fa1,2,ZHANG Shuai1,2

(1.School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China;2.Center for Engineering and Scientific Computation, Zhejiang University, Hangzhou 310027, China)

The surrogate modes have been wildly implemented in aerodynamic shape optimization. It is a fundamental problem to analyze the characteristic of surrogate modes for aerodynamic shape optimization. Four different surrogate models including a quadratic response surface and three Kriging surrogates were constructed based on the values of lift-to-drag ratio and volumetric efficiency responses. The structure of the various surrogate modes was investigated and the main differences were addressed. The Kriging surrogate based on exponential correlations produced a relatively better overall performance. The multi-object evolutionary genetic algorithm was applied to obtain the optimum design based on the surrogates for these two objectives. The two Kriging surrogate based on Exponential correlations and Gaussian correlations were selected and used in the optimization process. The results shows that the Kriging surrogate based on exponential correlations produces a better result for peak value.

surrogate mode;aerodynamic shape optimization;waverider configuration;aircraft conceptual design;multi-objective optimization

2015-01-21;

:2015-03-07。

國家自然科學基金(11102184)。

季廷煒(1982—),男,助理研究員,主要從事飛行器優(yōu)化設(shè)計方法研究。E-mail:zjjtw@zju.edu.cn

張帥,男, 副教授。E-mial: shuaizhang@zju.edu.cn

V211.5

A

1006-2793(2015)04-0451-07

10.7673/j.issn.1006-2793.2015.04.001

猜你喜歡
優(yōu)化模型設(shè)計
一半模型
超限高層建筑結(jié)構(gòu)設(shè)計與優(yōu)化思考
民用建筑防煙排煙設(shè)計優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
瞞天過海——仿生設(shè)計萌到家
設(shè)計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設(shè)計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
主站蜘蛛池模板: 亚洲第一视频网站| 亚洲成人一区二区三区| 精品無碼一區在線觀看 | 国产高颜值露脸在线观看| 91在线一9|永久视频在线| 亚洲精品无码人妻无码| 久久人妻xunleige无码| 黄色网站不卡无码| 伊人无码视屏| 狠狠v日韩v欧美v| 日本三区视频| 免费大黄网站在线观看| 成人无码一区二区三区视频在线观看| 香蕉久久国产精品免| 久久精品国产一区二区小说| 不卡视频国产| 亚洲欧美另类色图| 欧美啪啪一区| 国产喷水视频| 日本亚洲国产一区二区三区| 亚洲女同欧美在线| 40岁成熟女人牲交片免费| 日韩无码黄色网站| 国产欧美在线观看视频| 91久久偷偷做嫩草影院电| 国产欧美日韩精品综合在线| 久久精品一品道久久精品| 亚洲久悠悠色悠在线播放| 又大又硬又爽免费视频| 中文字幕2区| 国产美女精品人人做人人爽| 老司机精品99在线播放| 中文字幕人妻av一区二区| 伊人激情综合网| 国产又粗又猛又爽视频| 国产迷奸在线看| 亚洲娇小与黑人巨大交| 911亚洲精品| 亚洲av无码成人专区| 亚洲一区二区三区香蕉| 国模粉嫩小泬视频在线观看| 国产精品亚洲欧美日韩久久| 亚洲欧美成人网| 毛片网站在线播放| 亚洲自偷自拍另类小说| 亚洲最大综合网| 国产欧美日韩专区发布| 91免费精品国偷自产在线在线| 日本一区中文字幕最新在线| 99视频免费观看| 国产丝袜无码一区二区视频| 搞黄网站免费观看| 一级毛片基地| 色欲不卡无码一区二区| 亚洲欧美日韩久久精品| 99精品国产电影| 国产内射一区亚洲| 性激烈欧美三级在线播放| 色老头综合网| 91啦中文字幕| 又爽又大又黄a级毛片在线视频| 久久久亚洲色| 国产精品手机视频一区二区| 久久久久亚洲精品成人网| 久久国产精品77777| 久久天天躁狠狠躁夜夜躁| 亚洲国产精品美女| 一级毛片在线免费看| 久久国产亚洲欧美日韩精品| 日韩无码真实干出血视频| 97在线免费视频| 黄片一区二区三区| 国产精品播放| 网友自拍视频精品区| 日本不卡在线播放| 小说区 亚洲 自拍 另类| 日韩欧美国产三级| 精品成人一区二区| 国产乱子伦精品视频| 久久大香伊蕉在人线观看热2| 国产手机在线ΑⅤ片无码观看| 欧美亚洲国产日韩电影在线|