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

高效低噪的二維翼型優化設計

2017-03-09 10:14:21屈轉利李國才余培汛
振動與沖擊 2017年4期
關鍵詞:變形優化設計

李 鑫, 屈轉利, 李 耿, 劉 雨, 李國才, 余培汛

(1. 中國航天科技集團公司第四研究院第四十一所,西安 710025; 2. 西北工業大學 航空學院,西安 710072)

高效低噪的二維翼型優化設計

李 鑫1, 屈轉利1, 李 耿1, 劉 雨1, 李國才1, 余培汛2

(1. 中國航天科技集團公司第四研究院第四十一所,西安 710025; 2. 西北工業大學 航空學院,西安 710072)

針對傳統的氣動優化設計未考慮氣動噪聲影響的局限性,開展了基于噪聲預測模型的氣動優化設計方法在二維翼型中的應用研究。建立了由幾何外形參數化方法、徑向基函數(Radial Basis Function,RBF)動網格技術、改進粒子群優化算法、氣動分析方法、氣動噪聲預測方法等五大模塊構成的優化設計系統,且各模塊均采用標模算例進行驗證。通過對二維SC(2)-0714超臨界翼型進行了單點多目標優化設計。通過對比翼型幾何形狀、壓力系數分布以及在不同迎角下的氣動力系數曲線與總聲壓級的關系可得,翼型頭部半徑、厚度影響其頭部壓力峰值、壓力恢復、逆壓梯度等特性,從而影響升阻比和總聲壓級,逆壓梯度越小,翼型的總聲壓級越小。優化結果表明,在設計狀態下顯著提高了升阻比、降低氣動噪聲,考慮氣動噪聲的二維翼型優化設計系統可在實際的工程設計中進行應用。

外形參數化;RBF動網格;改進粒子群優化算法;氣動噪聲;翼型;優化設計

隨著航空業的發展,如何設計出更安靜的飛機是目前航空界共同關注的一個焦點。歐洲的ARACE[1](Advisory Council for Aeronautics Research)組織建議民用航空工業能夠在2020年前將乘客的人均噪聲污染降低到現在50%。由美國NASA[2](National Aeronautics and Space Administration)牽頭,FAA(Federal Aviation Administration)、飛機和發動機制造商攜手合作的“安靜技術驗證”項目,其目的就是將飛機的可感知噪聲水平在10年后降低50%,到未來25年后降低75%。由此可見,深入研究飛機噪聲的成因和傳播特性,進而提出改進方案以降低飛機的噪聲強度,具有重要的現實意義。如何在氣動外形優化設計中考慮氣動噪聲是一個興起的研究熱點。

另一方面,隨著計算流體力學技術和優化理論的不斷發展和完善,基于CFD(Computational Fluid Dynamics)技術的優化設計方法也逐漸成熟。該方法可以根據對飛行器性能的要求,直接以需要的氣動力、力矩等為設計目標,以一定的流體運動方程建立目標函數與飛行器幾何外形之間的關系,采用一定的優化算法,獲取目標函數的最優解。氣動優化設計方法最大的優點是使用靈活,受人為因素影響相對較小,適用于比較復雜的優化設計問題,在飛行器設計中廣泛應用[3]。國內外主要研究工作以考慮氣動力為設計目標,主要有:PALACIOS等[4]采用FFD方法對某超音速客機構型進行了氣動優化設計,顯著減小了激波阻力;黃江濤等[5]考慮了融合式翼梢小翼的穩健型優化設計;白俊強等[6]采用基于直接操作的FFD技術對RAE2822翼型進行了氣動減阻設計。然而,耦合氣動噪聲、氣動力的優化設計國內外研究工作比較少,主要有:POUANGUE等[7]研究了二維增升裝置考慮寬頻噪聲的優化設計;BIZZARRINI等[8]基于NAFNOISE求解器研究了高效、低噪的二維翼型優化設計 。

本文首先闡述了自由變形參數化、RBF動網格、改進粒子群優化算法、氣動力分析方法、噪聲預測模型等幾大部分的理論知識,并通過算例的測試對其進行了可靠性驗證。基于上述幾大模塊,建立了一套基于噪聲預測模型的翼型優化設計系統。最后,通過單點多目標約束對SC0014翼型進行優化設計,結果表明本文所采用的方法可為考慮噪聲的氣動優化設計提供了新的思路。

1 優化設計方法

1.1 自由變形參數化方法

自由變形技術(Free-Form Deformation,FFD)是一種氣動外形參數化技術,其優點在于穩健靈活,不需要對初始外形進行擬合,可以適用于空間任意幾何屬性的描述,且保持外形的幾何特性,如連續性、光滑性等。FFD方法是將研究對象(即待變形的翼型)嵌入FFD控制體中,通過移動各FFD控制點的位置來實現幾何外形的變形,因此在應用FFD方法進行幾何外形參數化的過程中,設計參數就是FFD控制點的位移。

在1986年首次提出FFD變形方法的論文中[9],通過Bernstein基函數來定義FFD控制體包圍的超曲空間幾何映射關系。首先在待變形幾何體周圍建立FFD控制體,然后計算待變形的幾何外形上每一個點在FFD 控制體中的局部坐標,局部坐標由式(1)確定

x(s,t,u)=

(1)

可以通過式(1)求出每個需要進行參數化的幾何外形上每個點的局部坐標。局部坐標對于給定的FFD控制點分布只需要計算一次,而不需要在每次幾何變形的過程中重復計算。在FFD控制點位置移動之后,FFD控制體中的待變形的幾何外形上任一局部坐標為 (s,t,u)的點x(s,t,u)的位移Δx由式(2)確定

Δx(s,t,u)=

(2)

(3)

變形后幾何外形上每一點的位置為

x′(s,t,u)=x(s,t,u)+Δx(s,t,u)

(4)

采用FFD方法進行翼型幾何外形參數化的過程如圖1所示,首先在初始翼型周圍建立FFD控制點(以黑色頂點表示),求解翼型每個幾何點在FFD控制點形成的控制體中的局部坐標(s,t,u)。然后以各個FFD控制點的位移為參數,改變FFD控制點的位置,在圖1中,拉動控制點使之位置發生改變。根據變形后的FFD控制點位置由式(2)計算變形后的翼型幾何外形,變形后的翼型如圖1(b)所示。圖2為本文所建立的FFD方法在飛機其它部件中的應用。

(a) 翼型變形前 (b) 翼型變形后圖1 FFD方法在翼型變形的應用Fig.1 Application of FFD method in airfoil

1.2 RBF動網格方法

在飛行器的優化設計過程中,其幾何外形不斷為了保證氣動分析的可靠性,需使計算域中的空間網格隨幾何外形的變化而自動更新,這一需求孕育出了動網格技術。由于RBF動網格方法[10]除變形后網格質量高之外,還具有適應性強、執行快速等特點,本文采用RBF方法生成中間過程中的計算網格。

構建思想是在初始網格的基礎上,建立邊界網格與空間網格之間的插值矩陣,空間網格的變形通過邊界網格的位移量來進行計算。其中邊界網格與空間網格之間插值矩陣的形式如下

(5)

(6)

(7)

下面分別以RAE2822翼型、M6機翼的網格變形來檢查程序的可靠性,具體情況分別如圖2所示。由圖可見,RBF動網格方法的魯棒性高,同時網格運動前后拓撲基本保持不變,也進一步保證了優化過程中CFD計算結果的可信度。對于4萬單元的計算網格,在優化之初構造表面網格和空間網格的插值矩陣耗時約1 min左右,以后每次重新建立空間網格的耗時不超過0.1 s,可見RBF動網格方法完全適用于翼型、機翼的優化設計中。

圖2 二維及三維動網格示意圖Fig.2 Mesh deformation of 2D and 3D

1.3 改進粒子群優化算法

優化算法的尋優能力、收斂速度都極大影響著優化設計效率。本文采用了基于二階震蕩和自然選擇的改進粒子群算法(Modified Particle Swarm Optimization,MPSO)[11],使用粒子二階震蕩來提高種群多樣性,并引入自然選擇思想來避免收斂速度過慢。其速度更新公式為

vi,j(t+1)=wvi,j(t)+

c1r1[pi,j-(1+ξ1)xi,j(t)+ξ1xi,j(t-1)]+

c2r2[pg,j-(1+ξ2)xi,j(t)+ξ2xi,j(t-1)]

(8)

文中采用二維函數Camel函數對該優化算法的計算效率進行了測試,其中Camel函數的表達式如下,x,y的取值范圍為-100~100

f(x,y)=

(9)

測試結果如表1所示,改進粒子群優化算法相對于粒子群優化算法效率更高,迭代次數更少。

表1 優化算法效率測試情況

1.4 氣動噪聲預測模型及驗證

Ffowcs Williams和Hall的翼型/機翼噪聲輻射理論已得到實驗和數值模擬的證實。根據FWH方程得到機翼后緣單元聲源遠場輻射噪聲強度的近似表達式為

(10)

式中:a∞為自由來流聲速;ρ∞為來流密度;v0為特征速度尺度;f0為特征頻率;L為翼型后緣距離觀測點的距離。考慮后緣角β的影響的情況下,根據Howe理論[12],噪聲強度的具體表達式為

(11)

式中:β為機翼后緣角;D(θ,ψ)為指向性項,指向性項表達式為

(12)

式中:ψ為方位角;θ為極角;具體意義如圖3所示。

圖3 指向性示意圖(方位角為90°)Fig.3 Sketch of direction

考慮到湍流斯特勞哈爾關系式為

(13)

式中,l0為湍流特征尺度。將式(13)代入式(11),則噪聲強度的表達式為

(14)

式中,v0和l0分別為沿zn方向的特征速度和特征尺度,其中zn為翼型上、下翼面的后緣垂直方向zn_up和zn_down。v0和l0由式(15)所求得,其中ω為最大湍動能TKE處的湍流耗散率,TKE和ω可從求解湍流模型為k-ω或k-ε方程的雷諾平均方程獲得,文中采用了湍流模型為k-ωSST方程的雷諾平均方程。

(15)

I可分為上下翼面Iup和Idown兩部分噪聲強度,I雖然不能準確表示噪聲的強弱,但它是一個相對準確值。可通過式(16)將I轉換成總聲壓級表達式

(16)

BROOKS等[13]通過噪聲風洞實驗,測量不同弦長、不同實驗狀態下的NACA0012翼型氣動噪聲,建立了一套翼型噪聲數據庫,并依此確立了一套半經驗翼型/機翼噪聲預測方法。本文選用BROOKS等用于試驗分析的二維NACA0012翼型,計算狀態與其試驗狀態相同,具體如圖4所示,并將計算結果與BROOKS等通過經驗公式計算的結果及ANOPP軟件計算的計算結果進行對比。定義噪聲監測點位于機翼后緣下方1.22 m,極角θ和方位角ψ均為90°。圖5為這3種方法的結果對比圖,圖5中NM及OASPL的定義為

NM(i)=10(0.1(NMi-NM1))

OASPL(i)=10(0.1(OASPLi-OASPL1))

(17)

式中,NMi、OASPLi分別指的是第i個工況情況下本文所采用的噪聲模型計算的結果、BROOKS等和ANOPP軟件計算的結果。由圖5可看出本文采用的噪聲模型所計算的結果與BROOKS等及ANOPP軟件計算結果相差很小,驗證了本文采用的噪聲預測模型具有一定的可靠性。

圖4 計算模型計算狀態圖Fig.4 State of computing models

圖5 計算結果對比圖Fig.5 Comparison of calculation results

由聲學理論可知,后緣(Trailing Edge,TE)噪聲與來流速度的五次方相關,而這個影響在本文采用的噪聲模型中也有所體現,因為特征速度u0將隨來流速度的改變而改變。另外模型中還包括了一些顯而易見的參數,例如指向性參數θ和ψ,后緣到觀測點距離H,后緣角β等。除了速度和這些參數外,在翼型外形優化中該模型還考慮翼型的幾何參數對噪聲計算的影響,主要體現在幾何參數變化帶來湍動能和湍流耗散率的變化。

2 基于噪聲預測模型的優化設計系統

基于第1.1節~第1.4節的理論內容,本節建立了基于噪聲預測模型的優化設計系統,分別由氣動外形參數化FFD模塊,網格自動生成模塊、氣動力分析模塊、氣動噪聲分析模塊、優化搜索模塊組成,其具體流程圖如圖6所示。

圖 6 基于噪聲預測模型的優化設計流程圖Fig.6 Optimization flowchart with noise prediction model

3 算例分析—單點多目標優化設計

在大型客機機翼設計過程中,翼型常采用超臨界翼型,為此本文選用了超臨界翼型SC(2)-0714為初始翼型。其優化設計狀態為:來流馬赫數0.2,單位尺度雷諾數4.61×106,迎角為8°,平均氣動弦長9.54 m,監測點位于后緣下方120 m處,極角和方位角為90°。優化設計目標為在起飛狀態下通過優化翼型提高其升阻比,降低噪聲,約束條件是翼型的最大厚度范圍保證13.5%~14.5%,定攻角8°。優化設計問題可表示為

在優化設計中,設置每代種群規模100個個體,進化代數60代后結果基本收斂,Pareto前沿線如圖7所示,其中縱坐標為總聲壓級差量(ΔNM=NMNow-NMoriginal,其中NMNow為優化歷程中當前翼型所計算得到的總聲壓級,NMoriginal為初始翼型的總聲壓級),橫坐標為升阻比。

圖7 Pareto前沿曲線Fig.7 Pareto frontier line

圖8 優化前后的翼型幾何對比Fig.8 Comparison of shape of original airfoil and optimized airfoil

優化前后的翼型幾何形狀、計算結果對比分別如圖8和表2所示,其中Original為初始翼型,K-max為升阻比最大的翼型,NM-min為總聲壓級最小的翼型,K-NM為介于兩優化目標間的一個翼型。對比翼型的幾何外形可以看出:相對于Original翼型,K-max翼型的彎度增大、上翼面頭部半徑增大、下翼面頭部半徑減小,而NM-min翼型其趨勢與K-max翼型相反。對于K-NM翼型來說,其翼型彎度、上翼面的頭部半徑相對于Original翼型增大,下翼面的頭部半徑減小,其趨勢與K-max翼型相似,但其厚度比K-max翼型要小。而由表2的計算結果可以看出,優化翼型K-max比初始翼型升阻比大1.246,總聲壓級降低0.104 dB;優化翼型NM-min相對于初始翼型來說,升阻比降低了0.58,總聲壓級降低0.679 dB;優化翼型K-NM相對于初始翼型來說,升阻比增大了1.06,總聲壓級降低0.51 dB。

表2 優化前后的翼型計算結果對比

圖9 聲壓級隨攻角的變化曲線Fig.9 Curves of SPL on different α

圖9為優化前后翼型在不同迎角下的總聲壓級、上下翼面所貢獻的聲壓級變化曲線圖,圖10為優化前后翼型升力系數隨迎角的變化曲線。結合圖9、圖10可以看出,當迎角為0°時,上翼面聲源對接收點處所貢獻的總聲壓級比下翼面要高, 主要是由于優化前后翼

型均具有一定的正彎度,雖然上下翼面均未出現分離,但由于彎度存在,上翼面的氣流順壓梯度相對于上翼面更弱,氣流脈動更強烈;隨著迎角的增大,升力系數升高,上翼面逆壓梯度增強,漸漸出現分離,導致上翼面聲源強度增強,對接收點的影響越來越大;隨著迎角的增大,升力系數增大,下翼面順壓梯度增強,抑制分離的出現,其聲源強度出現下幅度下降,相對于上表面其對接收點總聲壓級的貢獻變小;總聲壓級的變化主要來源于上翼面的聲壓級的貢獻,隨著迎角的增大,下翼面聲壓級的貢獻比重越來越小其大小。當出現大分離(α>12°)時,上翼面聲壓級急速增大。

圖11為優化前后翼型升阻比隨升力系數的變化曲線。從圖中可明顯的看出:優化后的翼型K-NM和K-max最大升阻比大于初始翼型。圖12為不同迎角下優化翼型K-NM與初始翼型Original的升阻比與總聲壓級之間的關系,黑色實線連著的兩個點表示同一迎角下的計算結果,從圖中可看出同一迎角下,K-NM的總聲壓級即小于初始翼型Original,且其升阻比也更大。通過對比分析可知該優化設計系統可進行高效低噪的二維翼型設計。

圖10 升力系數隨迎角的變化曲線Fig.10Liftcoefficientoverdifferentα圖11 升阻比隨迎角的變化曲線Fig.11Lift?dragratiooverdifferentα圖12 不同迎角下總聲壓級隨升阻比的變化Fig.12TotalSPLoverlift?dragratioondifferentα

為了進一步分析翼型的哪些因素對噪聲有影響,下面將從翼型壓力分布的角度分析與聲壓級的大小相關的因素。

文中選取了迎角8°時優化前后翼型的壓力系數分布對比曲線圖,具體如圖13所示。對比這幾組翼型可知:NM-min翼型上表面從40%位置處,其逆壓梯度最小,其次是K-NM、K-max、Original,這主要和翼型的彎度有關。而在這幾組翼型中,聲壓級的大小順序也是如此,文中猜測聲壓級的大小與翼型的彎度成正比關系,其主要是由于翼型的彎度越小,壓力恢復時,其出現的逆壓梯度越小,以致湍流的脈動越小,從而導致聲壓級較小;Original的負壓峰值更高,其次是NM-min、K-max、K-NM,這主要和翼型的頭部半徑大小相關。

4 結 論

通過采用所搭建考慮氣動噪聲的二維翼型優化設計系統,對比分析優化前后的翼型幾何形狀、壓力系數分布等特點,可得出以下結論:

(1)通過采用二維翼型的優化設計系統,對超臨界翼型SC(2)-0714進行考慮噪聲的外形設計,驗證了該系統可對比分析出不同計算狀態、不同幾何形狀、不同迎角下的翼型總聲壓級。

圖13 翼型壓力系數分布對比Fig.13 Comparison of the pressure coefficient distribution

(2)優化算例表明,應用該方法進行翼型的單點多目標優化設計,相比較于初始翼型,在不同迎角下,優化翼型K-NM在一定程度上降低氣動噪聲,提高升阻比。

(3)通過分析優化前后的翼型壓力系數分布可知:當翼型的彎度較小時,其對應的總聲壓級最小,總聲壓級的大小與翼型的逆壓梯度密切相關。

[ 1 ] HERKES B. The quiet technology demonstrator 2 flight test[C]//The Aviation Noise & Air Quality Symposium, 2006.

[ 2 ] DOBRZYNSKI W. Almost 40 years of airframe noise research-what did we achieve[J].Journal of Aircraft, 2010,47(2):353-367.

[ 3 ] 白俊強,劉南,邱亞松,等. 基于RBF動網格方法和改進粒子群優化算法的多段翼型優化[J]. 航空學報, 2013,34(12):2701-2715. BAI Junqiang, LIU Nan,QIU Yasong, et al. Optimization of multi-foil based on RBF mesh deformation method and modified particle swarm optimization algorithm[J]. Acta Aeronautica et Astronautica Sinica, 2013,34(12):2701-2715.

[ 4 ] PALACIOS F, ALONSO J J, COLONO M. Adjoin-based method for supersonic aircraft design using equivalent area distribution[C]//San Diego:50th AIAA Aerospace Sciences Meeting and Exhibit, 2012.

[ 5 ] 黃江濤,高正紅,白俊強,等. 基于任意空間屬性FFD技術的融合式翼梢小翼穩健型氣動優化設計[J]. 航空學報, 2013,34(1):37-45. HUANG Jiangtao, GAO Zhenghong, BAI Junqiang, et al. Study of robust winglet design based on arbitrary space shape FFD technique[J]. Acta Aeronautica et Astronautica Sinica, 2013,34(1):37-45.

[ 6 ] 白俊強,陳頌,華俊,等. 基于直接操作FFD技術的翼型氣動優化設計[J]. 航空計算技術,2013,43(1):40-43. BAI Junqiang, CHEN Song, HUA Jun, et al.Application of direct manipulated FFD technique in airfoil aerodynamic optimization[J]. Aeronautical Computing Technique,2013,43(1):40-43.

[ 7 ] POUANGUE A F, MNASRI C. Parameterization and optimization of broadband noise for high-lift devices[C]//Berlin:19th AIAA/CEAS Aeroacoustics Conference,2013.

[ 8 ] BIZZARRINI N, GRASSO F, COIRO D P. Numerical optimization for high efficiency, low noise airfoils[R]. AIAA Paper, AIAA-2011-3187.

[ 9 ] SEDERBERG T W, PARRY S R. Free-form deform of solid geometric models[J]. Computer Graphics,1986,20(4):151-160.

[10] ALLEN C B, RENDALL T C S. Unified Approach to CFD-CSD interpolation and mesh using radial basis functions[R]. AIAA Paper, AIAA 2007-3804.

[11] LI L, NIU B. The particle swarm optimization algorithm[M]. Beijing: Metallurgical Industry Press, 2009:25-29.

[12] HOWE M S. A review of the theory of trailing edge noise[J]. Journal of Sound and Vibration, 1978,61(3):437-465.

[13] BROOKS T F, POPE D S, MARCOLINI M A. Airfoil self-noise and prediction[M].Washington :NASA Reference Publication,1989.

A numerical optimization for high efficiency and low noise airfoils

LIXin1,QUZhuanli1,LIGeng1,LIUYu1,LIGuocai1,YUPeixun2

(1. The 41st Institute of the Fourth Academy of CASC, Xi’an 710025, China; 2. School of Aeronautics, Northwestern Polytechnic University, Xi’an 710072, China)

In view of the traditional aerodynamic optimization design without considering the effects of aerodynamic noise, the aerodynamic optimization design method based on noise prediction model was studied in the application of two-dimensional airfoil. The free deformation parametric method, radial basis function (RBF) mesh deformation technology, improved particle swarm optimization algorithm, pneumatic analysis method, and the aerodynamic noise prediction method for the optimization design, the five modules were established and formed a system. Finally, a single-point multi-objective optimization design of two-dimensional supercritical airfoil SC(2)-0714 was carried out. The effect of different airfoils geometry, pressure coefficient distribution, as well as the relationship between the aerodynamic coefficients and overall sound pressure level under different angles of attack were studied. The analysis showed that airfoil geometry could affect its head peak pressure, pressure recovery, adverse pressure gradient and other characteristics, which could change the lift-drag ratio and overall sound pressure level (SPL). And the adverse pressure gradient was closely related to overall SPL.Optimization results show that under the design condition the design system can significantly improve the lift-to-drag ratio, reduce the aerodynamic noise. It thus can be applied in practical engineering design.

parameterization method; radial basis function mesh deformation technology; particle swarm optimization; aeroacoustics; airfoil; optimization design

2015-11-03 修改稿收到日期:2016-01-24

李鑫 男,博士,工程師,1983年3月生

V211.3; V211.4

A

10.13465/j.cnki.jvs.2017.04.011

猜你喜歡
變形優化設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
“我”的變形計
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
例談拼圖與整式變形
主站蜘蛛池模板: 亚洲最大福利网站| 午夜精品一区二区蜜桃| 欧美精品影院| 91色国产在线| 久久狠狠色噜噜狠狠狠狠97视色| 亚洲免费成人网| v天堂中文在线| www亚洲天堂| 五月婷婷中文字幕| 茄子视频毛片免费观看| 国产精品自拍合集| 欧美日韩国产系列在线观看| 久久无码高潮喷水| 精品福利视频网| 99成人在线观看| 亚洲精品不卡午夜精品| 日日摸夜夜爽无码| 亚洲人人视频| 十八禁美女裸体网站| 狠狠色狠狠综合久久| 亚洲伊人天堂| 欧美日韩精品一区二区在线线| 国产欧美日韩专区发布| 国产一区二区三区日韩精品| 亚洲最大福利视频网| 亚洲综合激情另类专区| 天天综合网站| 欧美精品伊人久久| 亚洲成人黄色在线观看| 亚洲国产精品国自产拍A| 极品国产一区二区三区| 国产精品性| 亚洲精品中文字幕无乱码| 亚洲六月丁香六月婷婷蜜芽| www.亚洲国产| 日本高清免费一本在线观看 | 国产永久在线视频| 亚洲一级毛片免费观看| 久久精品视频亚洲| h网站在线播放| 国产97视频在线| 无码免费的亚洲视频| 欧洲亚洲一区| 日韩在线第三页| 亚洲天堂视频在线免费观看| 久久精品国产亚洲麻豆| 国产又爽又黄无遮挡免费观看| 真实国产乱子伦视频| 久久中文字幕不卡一二区| 中文字幕在线永久在线视频2020| 日本精品αv中文字幕| 素人激情视频福利| 久久91精品牛牛| 亚洲国模精品一区| 国产成人无码久久久久毛片| 丝袜高跟美脚国产1区| 亚洲欧美人成人让影院| 免费观看男人免费桶女人视频| 国产不卡在线看| 伊人精品成人久久综合| 国产农村1级毛片| 欧美成人看片一区二区三区| 青青久在线视频免费观看| 4虎影视国产在线观看精品| 国产激情在线视频| 亚洲第一黄色网| 国产精品漂亮美女在线观看| 亚洲va欧美va国产综合下载| 欧洲欧美人成免费全部视频| 欧美性精品| 婷婷丁香色| 看你懂的巨臀中文字幕一区二区| 无码国产偷倩在线播放老年人| 国产三级国产精品国产普男人| 亚洲天堂免费在线视频| 精品国产香蕉在线播出| 欧美一级高清视频在线播放| 国产制服丝袜91在线| 99re免费视频| 另类欧美日韩| 久久天天躁狠狠躁夜夜2020一| 亚洲a免费|