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

新概念融合升力體氣動布局設計優化方法研究

2015-06-26 15:48:32李治宇唐志共楊彥廣袁先旭唐小偉
空氣動力學學報 2015年1期
關鍵詞:參數化優化設計

李治宇,唐志共,楊彥廣,袁先旭,唐小偉

新概念融合升力體氣動布局設計優化方法研究

李治宇,唐志共,楊彥廣,袁先旭,唐小偉

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

參考國內外高升阻比飛行器氣動布局設計經驗,針對進出空間飛行器的氣動特性要求,開展跨速域高升阻比融合升力體氣動布局(BLB)研究以適應進出空間飛行器的各種要求,在傳統的翼/身外形的氣動效率與純升力體高容量效率之間尋求平衡。研究表明通過構建融合升力體數模,研究氣動外形的系統參數化描述方法,選擇設計變量及變化范圍,研究優化算法,建立融合升力體氣動布局設計及優化工具,開展融合升力體氣動外形優化設計是一種值得深入探討的研究方法。本文主要通過優化平臺集成數模參數化程序、網格自動化及基于Euler方程的快速流場求解程序進行優化設計并對優化結果進行分析計算,發展了一種快速有效的氣動布局優化設計方法,設計了初步滿足設計要求的新型高升阻比融合升力體氣動布局。設計的新布局能為再入飛行器氣動布局設計提供參考,所發展的優化設計方法計算速度快,成本低,可以為走向工程實用化的復雜外形氣動布局優化設計打下技術基礎。關鍵詞:融合升力體;參數化;升阻比;優化設計

0 引言

先進的氣動布局設計是各類高超聲速飛行器總體技術性能指標設計的關鍵所在,各國對此非常重視,并利用各種研究手段對其展開大量的研究。如美國航天飛機的布局設計,先后經歷了三個階段,設計了幾百個外形,最終選定現在的外形。目前各國的高超聲速飛行器氣動布局都以發展高升阻比、升力式再入復雜外形為特點,如乘波體、升力體、組合升力體和翼身組合體等。但升力體布局形式最主要的問題在于還不能較好地解決升阻比和防熱性能與操穩特性之間的矛盾。

HTV-2是美國研制的一種無動力高超聲速滑翔飛行器。在氣動布局上,采用了帶乘波體特征的高升阻比翼身融合體。按照設計要求,HTV-2將具有極高的升阻比,遠遠高于航天飛機,高升阻比使得彈道下降段航程更長,同時橫向機動范圍更大[1]。兩次失敗的飛行試驗表明,在氣動布局、操穩性能等方面的問題依然存在。

NASA的OSP(Orbital Space Plane)計劃指出為使再入飛行器具有較小的再入過載和較大的機動能力,要求飛行器高超聲速和亞聲速升阻比都要大;為了降低著陸速度,要求亞聲速升力系數較大等。該計劃提出的融合升力體氣動布局(Blended Lifting Body,BLB)試圖在傳統的翼/身外形的氣動效率與純升力體高容量效率之間尋求平衡,有望適應再入飛行器的各種要求[2-3]。

James S.Greathouse和Benjamin S.Kirk等人研究對比分析了返回艙外形、細長體、升力體以及有翼飛行器等各種外形的飛行器氣動力、熱性能,詳述了NASA選型參數指標。并運用牛頓流理論,應用簡單而準確的SNEWT程序進行了快速計算與比較分析。文獻中詳述了包括Overflow、CART3D、LAURA、Configuration Based Aerodynamics(CBAERO)等CFD代碼的特點和適用流動情況,并給出部分算例,在飛行器氣動外形優化設計的選型及計算方面給出參考[4-5]。

隨著我國航天技術的不斷發展,對各種飛行器的氣動性能提出了越來越高的要求,高超聲速飛行器的布局設計也越來越迫切。國內的多學科設計優化(MDO)理論研究尚處于起步階段,但MDO應用于設計研究的工作已經在各領域展開。國防科大的王振國等人在其專著中對多學科設計優化理論進行了詳細論述研究,舉例說明了多學科設計優化理論在飛行器設計中的應用[6]。西工大的江志國等人發展了自己的針對高超聲速飛行器氣動布局的優化軟件,總結出一套適用于高超聲速飛行器的性能指標計算模型、優化算法模型和優化思想,可應用于高超聲速飛行器的概念研究及初步設計[7]。車競等人將Pareto非劣解和遺傳算法結合起來,克服標準遺傳算法的缺點。為了提高優化效率,引入模擬退火算法進行局部搜索,形成了混合的模擬退火遺傳算法,并有效地用于高超聲速飛行器總體性能優化研究中[8]。唐偉、夏露等人也在高超聲速飛行器優化設計研究方面做出了大量的研究工作[9-10]。

由于國內尚未開展BLB概念的布局設計研究,本文擬參考美國BLB概念設計的經驗與方法及俄羅斯Clipper飛行器外形,通過對融合升力體布局形式和參數影響規律的研究,開展跨速域高升阻比融合升力體(BLB)氣動布局設計與優化,發展新型跨速域高升阻比再入飛行器。主要通過應用快速高效的優化設計平臺,計算分析特征尺寸對流場特性及飛行器氣動性能的影響,優化設計候選布局,并建立外形設計與優化的方法。這對未來發展高超聲速再入飛行器及相關氣動布局研究具有十分重要的意義。

1 研究方法

1.1 研究方法簡介

總體研究思路是基于優化平臺集成三維建模、網格自動化及氣動力快速計算模塊,形成優化設計方法。首先通過三維建模軟件作圖得到基本氣動布局外形并完成參數化建模,通過基于笛卡爾網格和Euler方程的流場快速求解方法計算氣動力系數,基于優化平臺選擇優化算法完成氣動布局優化設計,并針對優化得到的典型布局完成綜合氣動分析。

1.2 流場求解方法

在氣動布局優化設計中,快速準確地獲得氣動力特性非常重要。這是因為優化設計過程中的目標函數選取為升阻比等氣動特性,而這些氣動特性取決于飛行器的氣動布局和構形,并且與外型尺寸的關系無法用明確的函數表達出來。這就需要發展一種快速而準確的得到氣動力特性的方法,以滿足優化設計的精度和效率要求。

本文計算基于可壓縮Euler方程,采用有限體積法離散,計算網格使用直角笛卡爾網格,流動變量位于網格中心,時間離散采用Runge-Kutta法,通過時間推進得到穩態解,空間離散為迎風格式,選擇使用min mod限制器,格式具有TVD性質,采用多重網格法加速迭代過程。程序運行結束后輸出可讀氣動力系數文件[11]。

算例驗證顯示,該方法計算精度高,在高超聲速流動及亞跨聲速流動的數值計算中都有出色表現,可用于多種飛行器氣動力的計算分析。該方法的主要特點是計算速度快,效率較高[12]。數值計算中,網格量在130萬時,單機運行的計算時間可控制在10min之內。在優化設計時,在保證計算精度的前提下,適當減少網格量,可將單狀態運算時間控制在5min以內,能夠大幅縮短設計周期,降低計算成本。

1.3 參數化方法介紹

由于在優化設計過程中需要隨時生成更新的數模和網格,人工操作已經不能滿足計算要求。這就需要我們完成氣動外形的參數化,即發展一種能夠快速生成光滑的數模文件的方法。并且網格生成應該高度自動化,以便于提高計算效率,縮短優化周期。本文通過一個VB語言的腳本文件驅動外形特征變量的方法完成外形的參數化生成,并輸出parasolid格式的數模文件。在新外形生成時,驅動文件有錯誤參數的正確處理功能,如出現負值時,該文件并不會中斷程序,而是在輸出原始外形后繼續優化過程。所以,在出現錯誤尺寸時該文件不會導致優化過程中斷,也不會影響到優化結果,具有很好的魯棒性[11]。

2 優化算法驗證計算

2.1優化算法簡介

數值優化算法有快速高效、精度高的特點,但處理多峰值問題易陷于局部最優解,適用于單峰值優化問題。其基本思想就是迭代法,首先給出目標函數(設為f(x))的極值點的一個初始估計x(1)(稱為初始點),然后計算一系列的點x(2),x(3),x(4),……,希望該點列{x(k)}的極限值就是目標函數的極值點。公式表達為:

其中,S(k)為搜索方向,a(k)是一個正實數,稱為步長。不同數值算法的差別在于選取搜索方向和步長的方法不同,特別是選取搜索方向的方法。

本文選用直接數值優化技術中的修正可行方向算法MMFD(Modified Method of Feasible Directions)。該方法針對約束優化問題可以通過梯度搜索快速得到優化解,約束可以是等式或者不等式,在獲得最優解時可以以很高的精度滿足約束條件。MMFD遍尋當前設計點周圍區域,首先確定使得目標函數下降的搜索可行方向,沿著可行方向搜索;進行一維優化搜索,求得最優步長并確保下一步的優化點位于可行域內。在可行域內部搜索時,使用共軛梯度法確定搜索方向;在可行域邊界(即存在起作用約束時),搜索方

2.2 優化方法驗證計算

為驗證上述集成方法的有效性,通過對印度的RLV(Reusable Launch Vehicles)返回艙展開優化設計,并將優化結果與文獻數據對比,驗證其在雙目標優化設計中的應用[14]。RLV返回艙選擇球雙錐外形,外形輪廓及優化計算模型如下:

優化目標:min As,max Xcp

約束條件:V≥3m3

設計變量:0<Rn<1.0

選擇NCGA與MMFD相結合的優化算法進行優化[14],NCGA選擇群體大小30,進化代數25,交叉概率0.9,變異概率0.01。而文獻中同樣采用遺傳算法優化,其群體大小30,進化代數500,交叉概率0.9,變異概率0.001。優化結果對比如表1所示,外形如圖2。向應與起作用約束面成直角或者鈍角,即[gm(xk)]Tsk≤0。可行方向的確定方法是按梯度法搜索前進,這在不連續的優化空間里是難以實現的。

在優化過程中,為了避免陷入局部最優,在目標函數可能為多峰值時,需要采用全局優化算法。多目標優化設計中選用可以生成pareto解集的多目標優化算法——鄰域培植遺傳算法NCGA(Neighborhood Cultivation Genetic Algorithm)。在多目標優化設計中一般不存在多個目標同時達到最小的絕對最優解,不同目標之間的優化是相互沖突的,這就需要找到一個折衷的解集,在這個解集里不存在任何一個解可以在一個目標更加優化的同時而保證其他目標不變壞,這個解集即是pareto最優解,該解集中的解均為所設計問題的“非劣解”。領域培植遺傳算法是一種可以自動生成pareto最優解的遺傳算法,該算法可以通過調節遺傳操作參數控制優化過程,最終輸出pareto最優解和按給定目標權值得到的最優解[13]。

在多目標優化中得到的pareto最優解是一個由多個非劣解組成的解集,本文根據加權組合法從中選擇符合設計要求的最優解。加權組合法是將各自分目標函數fi(x)按照重要性分配權系數Wi(Weight)及比例因子Si(scalefactor),然后求和構成總的統一目標函數F(x),即:

圖1 印度RLV返回艙輪廓圖Fig.1Contour line of RLV reentry capsule

表1 印度RLV返回艙優化結果對比Table 1Optimization result for RLV reentry capsule

圖2 印度RLV優化pareto解集及外形比較Fig.2Pareto plot and the optimization contour of RLV reentry capsule

對比顯示,該方法能夠在滿足約束條件的情況下獲得比文獻中所列的優化結果更加優化的外形,即壓心位置更大,表面積和體積更小,達到了保證穩定性并且降低發射質量的優化目的,可以應用于飛行器氣動布局的多目標優化設計。

本文結果更優的原因在于遺傳算法的參數設置及數值優化算法的應用上。文獻[14]設置變異概率偏小,優化過程中遍尋設計空間,效率低,導致遺傳算法的優勢沒有發揮出來。另外,數值優化算法MMFD的應用也更進一步的優化了遺傳算法得到的最優解。

3 新型氣動布局優化設計

3.1 新型氣動布局設計及參數化

前期調研表明,俄羅斯Clipper飛船返回艙具有較高升阻比,在防熱性能及穩定性等方面也有較好表現[15-16]。參考Clipper外形,設計了新型融合升力體外形,外形及輪廓線如圖3所示。

該布局由三部分組成,分別是上部旋成體、中部弧-臺及下部弧面。

圖3 新布局外形及輪廓Fig.3New fashioned shape and contour

我們約定:圖中上部輪廓線各點標注為P1-P4,下部母線弧上的對應點,則相應標記為P11-P44,以示關聯,相對應的兩點橫坐標相等。以P1點為例,坐標表示為(XP1,YP1),依此類推。圓心坐標的下標以該圓的半徑表示,以下部母線大圓為例,其半徑為Rd,故圓心坐標表示為(XRd,YRd)。外形總長記為Length。原始外形的特征尺寸參數如表2所示。

以下表述各點及各特征尺寸間的幾何關系。

首先,Rn、L1、θ1、R1、L2、θ2、Rd、D為自變量。鑒于在后續計算中常用到各點坐標,故先推導圖中各點的坐標表示公式,如表3。

表2 輪廓線特征尺寸Table 2Characteristic dimension of contour line

表3 輪廓線特征點坐標Table 3Characteristic point coordinate of contour line

其中底部大圓圓心坐標(XRd,YRd)可以由下式推導得出:

以上所述,將外形輪廓線中所有特征尺寸表達出來,由8個自變量決定。

圖4 截面輪廓圖Fig.4Contour line of section

接下來介紹旋轉及放樣特征相關的幾何特征尺寸。如圖3所示,沿軸向分別在P2、P3、P4點取3個截面,編號為草圖2~草圖4。由于特征尺寸編號較多,為防止混淆,我們約定,尺寸下標只有一個數字的表示草圖編號,如D4代表草圖4所在的基準面到頂點的距離;尺寸下標兩個數字的,則第一個數字代表草圖編號,第二個數字為在本草圖內的編號,如R41代表草圖3中的標記為第1個的圓半徑;下標為數字和字母組合的,數字代表草圖編號,字母L代表左,R代表右(相對于左視圖方向),如R42R代表草圖4中右側第2個圓半徑。

其中,沿軸向的三個截面輪廓相似,由6個圓弧組成,如上圖4所示。這里以草圖4為例加以說明。草圖4中,自變量參數包括:θ41、R43、R44、θ44,應變量為R41、R42。其中R41=YP4,R42由幾個自變量共同決定。草圖3中,θ31、R33、R34與草圖4中相應特征尺寸相等,即草圖3中只增加了θ34一個自變量。同樣,草圖2中只增加了θ24一個自變量。這樣,由6個自變量參數可以完全描述沿軸向的三個截面。

綜上,該布局可以由14個自變量特征參數完全描述。

3.2 新型氣動布局優化計算

如前文所述,外形輪廓部分可由8個自變量參數描述,橫截面部分可由6個自變量參數完全描述。其中θ41、R43、R44為三個橫截面的共有參數,θ24、θ34、θ44為各截面自有參數。

作為初步優化設計,優化過程中保持基本輪廓形式不變,三個截面共有參數不變,即保持該布局頂部圓柱部分、過度圓弧半徑及底部大圓半徑不變,設計變量為θ24,θ34,θ44,在參數化外形合理并且保證能夠生成有效外形的范圍內約束設計變量變化范圍,并構建了基于原始外形的優化設計模型,優化目標為高超聲速高升阻比及進場亞聲速高升力系數最大化。優化模型可表示如下,其中的角度單位為弧度:

選用NCGA和MMFD相結合的混合優化算法,計算馬赫數分別取為:Mahyper=5.0,Masub=0.5,迎角均為10°。通過迭代優化計算,得到的優化外形如圖5所示。

圖5 優化結果示意圖Fig.5Schematic drawing of optimization result

3.3 優化氣動布局計算分析

本文完成了在馬赫數5.0和0.5情況下的多狀態氣動力計算分析,如圖6所示。

可以看到,該優化布局的高超聲速和亞聲速升阻比都比較高,迎角為0°時升阻比在1.0~1.3左右,并且隨著迎角的增加急速升高,在迎角8°~9°之間分別達到最大值,隨后隨著迎角的增加而降低。高超聲速最大升阻比約為2.14,亞聲速最大升阻比約2.95。

亞聲速升力系數較大,并且隨著迎角的增加而升高。迎角為0°時,亞聲速升力系數約為0.57,并且隨著迎角的增加而升高。在升阻比最大的8°迎角時,升力系數約為2.70。

由此可見,該優化布局在高超聲速再入及低速進場階段具有良好氣動特性,與Clipper原始外形[12]相比,同樣具有更好的升阻特性,是一種能夠初步滿足再入要求,值得深入研究的新型氣動布局。

在該例中,應用個人計算機運算完成了全部優化設計工作,節約了大量計算機時,降低了計算成本。

圖6 優化布局氣動力特性Fig.6Aerodynamic characteristics of optimization shape

4 結論

基于上述研究,通過新型融合升力體的氣動布局設計與優化,并對典型優化結果展開氣動特性計算分析,得到以下結論:

(1)設計的新型融合升力體氣動布局能夠初步滿足全速域高升阻比及進場亞聲速高升力系數的設計要求。優化設計結果較初始外形更“寬、扁”,在進一步的優化研究中應考慮容積率的影響,保證飛行器實用性。

(2)本文建立的氣動布局快速優化設計方法計算速度快,成本低,能夠應用于復雜外形飛行器的優化設計工作。

存在的不足之處主要表現在優化過程中沒有考慮氣動熱及操穩性能等方面的影響,變量約束還應再考慮容積率的影響。

下一步的工作將重點解決上述問題,進而考慮綜合氣動力、氣動熱、彈道及靜、動穩定性等方面,開展全面綜合的多約束多目標多變量的復雜外形優化設計工作。

[1]Henri D Fuhrmann.A blended lifting body aerodynamic design for the orbital space plane[R].AIAA 2003-3807.

[2]Douglas Nelson.Reentry performance analysis of an orbital space plane concept[J].AIAA 2003-7018.

[3]彭彪,張志峰,馬岑睿,等.國外高超聲速武器研究概述及展望[J].飛航導彈,2011,5:20-25.

[4]Benjamin M Andersen.Aerodynamic control on a lunar return capsule using trim-flaps[R].AIAA 2007-855.

[5]James S Greathouse,Benjamin S K,Randolph P L.Crew Exploration Vehicle(CEV)crew module shape selection analysis and CEV aeroscience project overview[R].AIAA 2007-603.

[6]王振國,陳小前,羅文彩,等.飛行器多學科設計優化理論和應用研究[M].北京:國防工業出版社,2006,04.

[7]江志國,唐碩,車競.高超聲速巡航飛行器氣動布局優化軟件設計[J].飛行力學,2008,26(1):52-55.

[8]車競,唐碩,何開鋒.高超聲速飛行器氣動布局總體性能優化設計研究[J].空氣動力學學報,2009,27(2):214-219.

[9]唐偉,高曉成,李為吉,等.雙橢圓截面再入飛行器的氣動計算及布局優化設計[J].空氣動力學學報,2004,22(2):171-174.

[10]夏露,高正紅,李天.飛行器外形多目標多學科綜合優化設計方法研究[J].空氣動力學學報,2003,23(1):275-280.

[11]李治宇.載人登月返回艙氣動布局概念設計與優化方法研究[D].[碩士學位論文].綿陽:中國空氣動力研究與發展中心,2011.06.

[12]李治宇,袁先旭,楊彥廣,等.基于笛卡爾網格的無粘氣動力計算[C]//第十五屆全國計算流體力學會議,2012.08.

[13]周明,孫樹棟.遺傳算法原理及其應用[M].北京:國防工業出版社,1999.06.

[14]Rajesh Kumar Arora.Aerodynamic shape optimization of a reentry capsule[R].AIAA 2003-5394.

[15]Pavel Vashchenkov,Andrey Krylov.Numerical analysis of high altitude aerodynamic of reentry vehicles[C]//13th International Space Planes and Hypersonic Systems and Technology.AIAA 2005-3409.

[16]Mikhail S I,et al.Numerical modeling of high altitude aerodynamic of reentry vehicles[C]//9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference.AIAA 2006-3801.

[17]Alex Van der Velden,et al.iSIGHT-FD,a Tool for multi-objective data analysis[R].AIAA 2008-5988.

Aerodynamic shape design and optimization method for a new blended lifting body

Li Zhiyu,Tang Zhigong,Yang Yanguang,Yuan Xianxu,Tang Xiaowei
(China Aerodynamics Research and Development Center,Mianyang621000,China)

Taking the references of the high lift-to-drag ratio vehicle aerodynamic shape design experiences,aimed at the aerodynamic performance requirements of the reenter vehicle,a blended lifting body aerodynamic shape with high lift-to-drag ratio is investigated to fit the various requirements of the reentry vehicle,and to quest a balance between the aerodynamic performance of the traditional wing/body shape and the high cubage efficiency of lifting body.The prophase research shows that it is an advisable method to develop the blended lifting body aerodynamic shape optimization method by building the blended lifting body model,researching the parameterization method for the blended lifting body aerodynamic shape,selecting the design variables and their range,researching the optimization algorithm,building the blended lifting body aerodynamic shape design and optimization tool.In this paper,the optimization process is proposed by integrating the geometry dimensions parameterization program,the automatic grid generation technique and the fast aerodynamic numerical calculation program based on Euler equations.A fast aerodynamic shape optimization method is developed,and the optimization result is obtained and analyzed.A new blended lifting body shape with high lift-to-drag ratio is designed.The new shape provides an useful reference for the reentry vehicle design,and the optimization method has obvious advantages in computing speed and cost,and lays a technique foundation for the practical complex aerodynamic shape optimization.

blended lifting body;parameterization;lift-to-drag ratio;optimization

V211.3

Adoi:10.7638/kqdlxxb-2013.0047

0258-1825(2015)01-0048-06

2013-04-22;

2013-06-10

李治宇(1985-),男,碩士,助理工程師,研究方向:氣動布局優化設計及高超聲速氣動力計算.E-mail:lizhiyu@mail.ustc.edu.cn

李治宇,唐志共,楊彥廣,等.新概念融合升力體氣動布局設計優化方法研究[J].空氣動力學學報,2015,33(1):48-53.

10.7638/kqdlxxb-2013.0047.Li Z Y,Tang Z G,Yang Y G,et al.Aerodynamic shape design and optimization method for a new blended lifting body[J].Acta Aerodynamica Sinica,2015,33(1):48-53.

猜你喜歡
參數化優化設計
從一道考研題談空間曲線積分的計算
Pro/E的三維往復壓縮機參數化模型庫的建立
一種懸架運動仿真快速建模方法研究
汽車科技(2016年6期)2016-12-19 20:32:56
汽車行李箱蓋鉸鏈機構的分析及優化
東林煤礦保護層開采卸壓瓦斯抽采優化設計
橋式起重機主梁結構分析和優化設計
對無線傳感器網絡MAC層協議優化的研究與設計
科技視界(2016年22期)2016-10-18 15:25:08
基于simulation的醫用升降椅參數化設計
科技視界(2016年21期)2016-10-17 17:27:09
簡述建筑結構設計中的優化策略
股骨頸骨折內固定方式優選方法研究
主站蜘蛛池模板: 久久久久亚洲AV成人网站软件| 亚洲妓女综合网995久久| 亚洲色图另类| 最近最新中文字幕在线第一页| 尤物精品国产福利网站| 99热国产这里只有精品无卡顿"| 亚洲成人黄色在线观看| 亚洲精品男人天堂| 五月婷婷欧美| 91精品小视频| 波多野结衣无码中文字幕在线观看一区二区 | 少妇精品网站| 日韩A∨精品日韩精品无码| 国产00高中生在线播放| 91九色国产porny| 九九热免费在线视频| 国国产a国产片免费麻豆| 18禁黄无遮挡免费动漫网站| 亚欧美国产综合| 成人午夜视频免费看欧美| 在线日韩日本国产亚洲| 永久成人无码激情视频免费| 激情六月丁香婷婷四房播| 久久久久国产精品熟女影院| 四虎成人精品| 精品国产自在在线在线观看| 精品色综合| 在线精品欧美日韩| 动漫精品啪啪一区二区三区| 午夜少妇精品视频小电影| 国产午夜不卡| 五月婷婷伊人网| 91麻豆精品国产高清在线| 免费一看一级毛片| 国产对白刺激真实精品91| 一级毛片高清| 亚洲精品国偷自产在线91正片| 污污网站在线观看| 亚洲欧美日韩另类| 久热中文字幕在线| 欧美色综合久久| 日韩精品一区二区三区免费| 天天摸天天操免费播放小视频| 亚洲精品免费网站| 欧美激情一区二区三区成人| 丁香五月婷婷激情基地| 国产日本一线在线观看免费| 国产不卡网| 欧美在线视频不卡| 日本免费a视频| 亚洲高清国产拍精品26u| 欧美不卡二区| 综合色在线| 另类欧美日韩| 亚洲天堂日韩在线| 国产在线八区| 一级毛片免费的| 国产尹人香蕉综合在线电影| 色噜噜久久| www.youjizz.com久久| 91外围女在线观看| 久久免费观看视频| 日韩久久精品无码aV| 无码高潮喷水在线观看| 国产国产人成免费视频77777| 重口调教一区二区视频| 国产成人精品18| 午夜日本永久乱码免费播放片| 色妺妺在线视频喷水| 亚洲欧美人成电影在线观看| 真人免费一级毛片一区二区| 有专无码视频| 免费大黄网站在线观看| 456亚洲人成高清在线| 六月婷婷综合| 亚洲精品国偷自产在线91正片| 国产成人做受免费视频 | 精品国产成人高清在线| 亚洲一区二区三区麻豆| 亚洲黄网在线| 一区二区三区毛片无码| 国产女人18水真多毛片18精品 |