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

水平風洞中開展飛機尾旋特性研究的理論分析

2018-07-23 09:15:14郭林亮祝明紅傅澔楊洪森鐘誠文
航空學報 2018年6期
關鍵詞:飛機模型

郭林亮,祝明紅,傅澔,楊洪森,鐘誠文

1. 西北工業大學 航空學院,西安 710072 2. 中國空氣動力研究與發展中心 低速空氣動力研究所,綿陽 621000

尾旋作為飛機最復雜、最危險的極限飛行狀態之一,以超臨界迎角、大側滑角、顯著的滾轉角速度及偏航角速度為主要特征。由于尾旋中飛機的操縱性顯著變壞甚至完全喪失,同時空間方位判斷及駕駛條件變得更為復雜,因而對飛行員及飛機的安全構成了極大威脅。因此,飛機的大迎角失速偏離、尾旋特性預測以及尾旋改出方法,都是飛機氣動設計和地面試驗中非常重要的研究課題。

現有飛機尾旋相關的研究手段主要有飛行仿真和試驗兩大類。其中飛行仿真手段側重于解析分析,主要包括非線性分叉分析、延拓算法研究、數值飛行動力學仿真及基于飛行模擬器的半物理仿真等,其所需的氣動力數據由風洞試驗獲得,主要有大迎角靜態試驗、動導數試驗、旋轉天平試驗3類[1-3]。近年來,NASA開展了縮比大型運輸機的氣動力建模、分叉分析仿真及尾旋試驗研究[4],歐盟第七框架項目SUPRA開展了大型飛機針對極限飛行條件下移動基模擬器的飛行員評價仿真研究[5-6]。理論研究中一些新的尾旋敏感性分析方法相繼提出,如SNB(Saddle-Node Bifurcation)準則[7]、非線性指數理論(Nonlinearity Index Theory)[8],這些方法為飛機尾旋預測提供了新的途徑,但需要相應試驗方法的驗證和支持。試驗手段主要包括立式風洞尾旋試驗、縮比模型大氣自由飛試驗以及原型機的尾旋試飛等[9-10]。立式風洞尾旋試驗主要研究飛機穩定尾旋特性及改出操縱方法,其環境可控、成本低廉、效率較高,但無法模擬尾旋進入階段;縮比模型的大氣自由飛試驗可進行尾旋進入、發展和改出全過程的研究,但系統復雜、耗費較高、效率有限;原型機試飛是為設計定型開展的驗證飛行,其結果真實可靠,一般僅在設計研制過程的最后階段開展。

近年來發展的風洞虛擬飛行技術是典型的針對非線性氣動問題的試驗手段,能夠較為逼真地模擬機動飛行過程,更好地揭示氣動/運動的耦合機理。該技術最早由美國人提出概念并由Magill等在喬治亞州42 in×43 in(1 in=25.4 mm)風洞進行了首次虛擬飛行的演示試驗[11-14],隨后其他國家也開展了類似的研究。英國Bristol大學設計了低速風洞3自由度和5自由度動態裝置,發現了Hawk模型大迎角下的極限環振蕩現象,并實現了該模型繞速度矢滾轉動作[15-19]。俄羅斯Sohi發展的腹撐3自由度裝置開展了典型戰斗機穩定尾旋的試驗研究,但沒有對系統動力學進行理論建模分析,機構中的運動曲桿和摩擦力矩對試驗結果的影響尚未涉及[20]。俄羅斯TsAGI提出了一種背撐3自由度裝置,采用魯棒控制方法對機翼搖滾問題進行了抑制研究[21]。中國空氣動力研究與發展中心在2.4 m跨聲速風洞開展了高機動導彈的虛擬飛行試驗[22],在?3.2 m低速風洞開展了針對固定翼飛機的虛擬飛行試驗,目前試驗達到的最大迎角接近30°[23]。

綜上分析,風洞虛擬飛行試驗技術具備多自由度、角度范圍大、模型自由轉動等特點,可用于飛機研制初期開展大迎角及失速尾旋特性研究。本文提出了一種大角度虛擬飛行試驗裝置,能夠開展大迎角失速/偏離/尾旋全過程的模擬研究,可成為立式風洞尾旋試驗的有力補充,可進一步完善針對尾旋問題的動態試驗技術研究體系。本文采用拉格朗日乘子法建立了該裝置的數學模型,結合已有風洞試驗數據開展仿真研究,并與立式風洞尾旋試驗結果進行對比分析,驗證了該技術的可行性和應用潛力。

1 試驗裝置

1.1 尾旋運動特征分析

圖1 穩定尾旋的受力分析Fig.1 Force analysis of steady spin

尾旋是飛機在超臨界迎角范圍出現繞其縱軸的自轉后,在氣動力、慣性力及重力的作用下,一方面繞其自身3個體軸旋轉,另一方面沿半徑很小的螺旋形軌跡做自發的下降運動,如圖1所示。在尾旋的初始階段,由于大迎角非對稱力矩或舵面產生的偏航力矩使得飛機開始旋轉;之后由于慣性耦合作用,飛機迎角進一步增大從而進入穩定尾旋狀態。穩定尾旋中,重力和阻力基本平衡,升力提供向心力以維持繞鉛垂軸的穩定旋轉;三軸力矩基本保持平衡。在多數情況下,尾旋半徑對飛機尾旋運動參數的影響較小,在一段時間內高度變化帶來空氣密度的影響也可忽略。因此飛機的尾旋運動可用一個動力學相似縮比模型繞一固定軸的3自由度旋轉來進行模擬。比如現有的旋轉天平試驗、立式風洞尾旋試驗等。

1.2 裝置設計

本文提出的3自由度裝置將在水平風洞中實現接近90°迎角的飛行動作模擬,可開展失速偏離、尾旋初始階段和穩定尾旋的研究,如圖2所示。具體的技術方案是:翼型垂直支桿通過下端安裝面固定在低速風洞下洞壁或平臺上。輕質曲桿一端與旋轉鉸聯結,繞平行于風洞中心軸線自由旋轉;另一端與十字萬向鉸固聯,萬向鉸置于模型內部用來支撐模型,通過旋轉鉸和萬向鉸的共同作用可實現模型在俯仰、偏航、滾轉3個方向較大范圍內的自由轉動。位于旋轉鉸前端的引電器用于模型內部的機載設備供電及采集機載設備的數據信號,與曲桿同步旋轉,如圖3所示。曲桿端部配重的高度可以手動調節,以保證曲桿的重心在曲桿的旋轉軸上,保持試驗旋轉運動過程中實現動平衡。這樣通過滾轉軸和俯仰/偏航二自由度轉臺相結合的方式實現了試驗模型物理意義明晰的3自由度運動,可實現尾旋進入、發展和改出全過程的模擬。

圖2 虛擬飛行裝置示意圖Fig.2 Schematic of virtual flight test rig

圖3 虛擬飛行裝置繞速度矢滾轉示意圖Fig.3 Sketch of velocity-vector roll in virtual flight test rig

圖4 虛擬飛行試驗裝置變形云圖Fig.4 Strain contours of virtual flight test rig

圖5 虛擬飛行試驗裝置模態振型Fig.5 Vibration mode of virtual flight test rig

圖4給出了虛擬飛行試驗裝置的變形云圖,施加載荷為后文試驗模型X、Y、Z3個方向上的最大氣動載荷,此時曲桿末端的最大位移約為3.5 mm,變形量基本可以接受,后續可通過優化設計進一步提高剛度。圖5給出了模態振型分析結果,其一階固有頻率f約為7.6 Hz,結合后文試驗模型的動態響應特性,應可避免發生共振,圖中色例描述的是歸一化的位移量。

2 動力學模型

2.1 數學模型建立

2.1.1 前提和假設

風洞試驗時,模型質心位于風洞試驗段中心保持不變,但模型的姿態可繞三軸轉動,因此,在動力學建模時可作如下假設:

1) 假設飛機模型及試驗裝置均為剛體,不存在變形情況。

2) 假定曲桿氣動力為零,曲桿的旋轉軸平行于風洞來流方向。

3) 模型的轉動中心為其重心,且模型和曲桿的重心在曲桿的旋轉軸上。

4) 試驗模型的慣量和旋轉速率較小,可忽略陀螺力矩的影響。

2.1.2 軸系及角度定義

為方便建立數學模型,本文使用如下坐標系及角度定義(見圖6):

地軸系OgXgYgZg,簡稱Sg,其OgZg軸沿鉛垂方向向下,OgXg軸在水平面內,并與風洞軸線平行,與來流速度V方向相反,OgYg軸與平面OgXgZg垂直,指向右。

模型體軸系OaXaYaZa,簡稱Sa,原點為模型重心,OaXa軸在飛行器對稱平面內,平行于機身軸線或機翼的平均氣動弦線,指向機頭;OaYa軸垂直于對稱平面,指向右;OaZa軸位于對稱面內,指向符合右手定則。

圖6 坐標系定義Fig.6 Definition of coordinate system

曲桿體軸系OrXrYrZr,簡稱Sr,原點位于曲桿重心,軸系定義與Sa類似。

2.1.3 基于絕對坐標方法的動力學方程

按照絕對坐標方法,定義本系統的絕對坐標為

同時定義基于自身體軸的角速度向量ωi=[piqiri]T(i=r,a)。并有如下關系:

(1)

(2)

根據拉格朗日乘子法,給出動力學方程的一般形式為[24]

(3)

式中:λ為拉格朗日乘子;A為與無約束動力學方程有關的系數矩陣,即

(4)

其中:Ja、Jr分別為飛機模型和曲桿的慣量矩陣;Da、Dr分別為飛機模型、曲桿體軸系角速率和地軸系姿態變化率的關系矩陣。

B為該方程組右側向量,即

(5)

(6)

2.1.4 無違約算法

多體系統動力學建模時,常會遇到約束條件違約的問題,特別是當模型進入振蕩過程中,違約現象比較明顯。通常情況下這是由于動力學模型中對于約束條件的選取往往只對位移、速度、加速度中的一個量進行約束,在數值計算中不可避免地會出現其他量違背約束條件的情況。

(7)

(8)

此時將速度約束表達式引入加速度約束表達式,得到

(9)

于是原動力學方程式(3)改寫為

(10)

以萬向鉸幾何約束為例,考核原動力學方程式(3)與無違約算法之間的優劣。圖7是在相同條件下兩種算法對萬向鉸兩軸夾角進行復算的結果對比。如圖所示,原動力學方程的結果顯示萬向鉸兩個旋轉軸夾角不斷擴大,與實際情況不符;采用無違約算法之后,幾何違約情況基本消除,偏差值極小,修正效果好。

圖7 違約修正影響Fig.7 Effect of constraint violation correction

2.2 氣動力模型

本文基于常規風洞試驗得到的基本氣動力、舵面效率、動導數及旋轉天平等數據建立了如下氣動力模型,各系數均在體軸系下描述。

(FA,MA)=f(α,β,ωa,V,δe,δr,δf,δa)

(11)

式中:α為迎角;β為側滑角;ωa為飛機的旋轉角速度;δe、δr、δf、δa分別為模型升降舵、方向舵、襟翼和副翼角度,各氣動系數按照基本量+增量的方式計算,即

ΔCiFO(po,qo,ro)+ΔCiRB(α,β,ωss)

i=l,m,n;j=a,e,r

(12)

式中:CiST為與迎角α、側滑角β、襟翼δf有關的靜態基本量;下標l、m、n分別表示滾轉力矩、俯仰力矩和偏航力矩;ΔCiδj(δj)為舵面δj帶來的增量;ΔCiFO為機體轉動引起的動導數增量;ΔCiRB為機體繞速度矢旋轉產生的增量;ωss為機體繞速度矢旋轉的速度,即

ωss=pbcosαcosβ+qbsinβ+rbsinαcosβ

(13)

其中:pb、qb、rb為體軸系下的三軸角速率;posc、qosc、rosc用于計算動導數的振蕩分量,其表達式為

(14)

2.3 摩擦力矩估算

常用的動態摩擦模型有Dahl模型、Bristle模型、Bliman-Sorine模型、LuGre模型和Leuven模型[25-26]。其中LuGre模型能精確地描述摩擦靜態和動態特性,因此選擇LuGre模型在穩態系統下的一種簡化形式對裝置中兩組鉸摩擦力矩進行計算,具體形式為

(15)

2.4 質心偏移影響

理想狀態下,模型質心與虛擬飛行支撐裝置萬向鉸的旋轉中心應該重合。但實際上,為了保證無風條件下系統的穩定性,模型設計時質心應略低于旋轉中心。另外,模型機身的高度有限,其Z向上的質心偏移較其余兩個方向也不易調整。

質心偏移對試驗系統帶來的影響主要有3個方面:模型慣量變化、重力引起的力矩變化、氣動力產生的力矩變化。對于模型慣量的變化,按照慣量平行移軸進行處理;而對于后兩者,將重力、氣動力在模型體軸系表述后,與偏移向量叉乘可得到力矩的變化量,具體為

MG=rcg×G

(16)

MFA=rcg×FA

(17)

式中:MG、MFA分別為重力和氣動力產生的力矩;rcg為質心偏移向量;G、FA分別為重力和氣動力。

3 仿真與尾旋試驗方法

3.1 仿真程序

仿真程序在MATLAB環境下編寫,采用Nelder-Mead算法在給定初始狀態下配平飛機;導出約束方程的Jacobian矩陣的解析形式,與無約束的動力學方程聯立構成動力學模型,然后采用無違約算法進行求解。圖8給出了仿真程序的相關模塊和主要流程。

圖8 仿真流程Fig.8 Flow chart of simulation

3.2 立式風洞尾旋試驗

圖9 立式風洞尾旋試驗場景Fig.9 Scene of spin test in vertical wind tunnel

為了驗證該方案的可行性,將仿真結果與立式風洞試驗結果進行對比分析。?5 m立式風洞是國內唯一一座研究飛機尾旋運動的特種風洞。尾旋試驗模型除與真實飛行器保持幾何相似外,還需要滿足動力學相似準則,主要是滿足弗勞德數相似[27-28]。立式風洞采用吊掛支持模型法進行飛機尾旋研究試驗。試驗中,飛機模型被兩根繩索懸掛在風洞試驗段中,如圖9所示。在試驗段氣流較低時,模型投手將模型按預定的姿態(模型迎角大于失速迎角)投入試驗段中心區域垂直上升的氣流中,同時給模型施加一個初始旋轉速度(約1~3圈/s),然后逐步增加試驗段風速到某個值使得模型平衡于風洞中心區域并進入尾旋狀態(平衡燈亮),尾旋運動測量系統開始記錄模型做尾旋運動時各運動參數隨時間的變化規律。模型尾旋運動參數由微型姿態航向參考系統(Attitude Head Reference System, AHRS)獲得,AHRS內部集成了基于MEMS(Micro ElectroMechanical Systems)技術的三軸陀螺儀、加速度計及磁強計,其內部包含了嵌入式的姿態數據解算單元,通過Kalman濾波的方法,給出模型姿態角、角速度、加速度等測量信息。

3.3 仿真與試驗條件

本文以某飛機為例,按照相同的進入舵偏、止旋舵偏、低頭舵偏開展研究,仿真中飛機模型的主要參數與立式風洞尾旋試驗的模型參數保持一致,具體見表1。

首先通過與立式風洞尾旋試驗結果、6自由度動力學方程仿真結果進行對比,驗證3自由度虛擬飛行試驗裝置數學模型的正確性和合理性。仿真與試驗的舵面控制信號均為:進入尾旋階段,1 s時升降舵為-5°,2.5 s時升降舵回零,副翼為30°;改出尾旋階段,18.4 s時使副翼偏轉到-30°,1 s后副翼回零,升降舵為30°。

在此基礎上,進一步分析了該裝置中曲桿慣量、摩擦力矩以及模型質心偏移等因素對結果的影響,給出設計機構和開展試驗的建議。此時的控制信號僅保留上述進入尾旋階段的信號。

表1 試驗模型主要參數Table 1 Major parameters of test model

4 結果分析

4.1 與立式風洞尾旋試驗結果對比

如圖10所示,圖中3組曲線分別是6自由度仿真計算結果、3自由度曲桿虛擬飛行裝置無摩擦仿真結果(圖中標注為“帶曲桿仿真”)和立式風洞尾旋試驗結果。由于立式風洞尾旋試驗不能模擬尾旋進入階段,因此,該曲線僅有穩定尾旋及改出階段。6自由度和3自由度虛擬飛行的結果總體一致性較好,主要區別在于尾旋進入階段,前者的俯仰振蕩較后者劇烈一些,這可能與3自由度虛擬飛行裝置的位移約束有關。

從迎角響應可以看出,尾旋穩定階段三者的平均迎角均在82°左右,振蕩周期接近,其中3自由度虛擬飛行的結果略大。3者的振蕩谷值比較接近,但兩組仿真數據曲線比試驗結果振幅小。采用同樣的止旋操作、低頭操作后,數據響應總體趨勢一致,3種狀態均能達到改出效果。同時注意到,尾旋試驗中迎角低于60°后不再下降,這與尾旋試驗時一旦改出尾旋狀態后模型迅速脫離流場區域以及上下吊索的拖拽有關。從側滑角響應看出,3者的平均側滑角均在7°左右,振蕩周期3者接近,兩組仿真結果的振幅同樣比試驗結果小,改出過程趨勢基本一致。

6自由度和3自由度的仿真結果與立式風洞試驗結果在迎角、側滑角振蕩幅值等方面的差異,一方面與仿真中使用的靜態測力、動導數及旋轉天平等風洞試驗數據的系統誤差、動態試驗數據使用方法有關,試驗誤差包括洞壁干擾、支架干擾、雷諾數效應等;另一方面與立式風洞尾旋試驗的誤差有關,包括數據精度、上下吊索干擾、模型可能存在豎向的重心偏移等因素。

從三軸角速度響應看出,仿真數據與試驗數據的均值接近,俯仰和滾轉振蕩幅值較小,改出階段響應趨勢一致。3組結果的繞速度矢滾轉速率比較接近,其中3自由度虛擬飛行的結果略小,改出過程變化趨勢基本一致。

此外,虛擬飛行試驗中,由于風速始終保持水平方向,模型繞速度矢滾轉時模型俯仰姿態呈正弦振蕩形式,但基本不影響氣動力。3者穩定尾旋狀態的下沉速度基本一致。

綜合以上對比分析認為,3自由度虛擬飛行試驗裝置可以模擬飛機尾旋進入、發展和改出的全過程;與立式風洞試驗結果對比表明,該技術可以捕獲尾旋發展和改出階段的主要特征參數,為風洞中研究大迎角失速/尾旋問題提供新的解決途徑。

圖10 仿真結果與垂直風洞試驗數據對比Fig.10 Comparison of results of simulation and vertical wind tunnel test

4.2 旋轉曲桿的影響

3自由度虛擬飛行試驗裝置中模型的滾轉運動通過曲桿和模型的同步旋轉實現,因此曲桿的存在對模型的運動會有一定的影響。由于曲桿只有一個轉動自由度,其繞自身X軸的慣量是影響模型尾旋運動的主要因素。

圖11給出了無摩擦狀態下不同曲桿慣量的仿真結果。本文中默認的曲桿慣量占飛機模型滾轉慣量的10%,以此為基準,分別減小和增大50%,得到3組不同慣量下的結果。從圖中可以發現,模型經過一段時間的大幅振蕩后才進入穩定尾旋狀態,而隨著曲桿慣量增大,進入階段的時間明顯增加,幅值有所增大,但振蕩周期基本不變。進入穩定尾旋后,迎角、側滑角振蕩的幅值隨著曲桿慣量的增加略有增長,振蕩均值基本不變,振蕩周期基本也不變;角速率方面受曲桿慣量的影響與之類似,振蕩幅值增加、均值和周期基本不變,繞速度矢滾轉速率基本不變。因此曲桿對穩定尾旋影響有限,但對進入階段有一定影響。

通過對迎角、側滑角穩定振蕩的部分進行快速傅里葉變換(FFT)并繪制頻譜(見圖12)后發現,模型迎角、側滑角振蕩存在兩個主要頻率,分別為1.1 Hz和1.35 Hz,其中后者與尾旋自由飛時模型振蕩頻率一致,因此前者應為支撐機構帶來的影響。此時兩個頻率相近的低頻信號互相疊加,引發拍現象,使得動態響應曲線的振蕩幅值忽大忽小。

圖11 不同桿慣量仿真結果Fig.11 Results of simulation with different inertia of rod

圖12 不同桿慣量仿真結果的幅頻圖Fig.12 Amplitude-frequency diagram of results of simulation with different inertia of rod

從頻譜圖亦可看出,隨曲桿慣量增長,3組數據的頻率基本不變。這是由于穩定尾旋發生在模型迎角82°附近,此時將曲桿慣量投影到模型體軸系后發現,曲桿慣量分解到模型偏航軸的部分占主量;而模型偏航慣量近似為滾轉慣量的8倍,經過測算,曲桿慣量的變化使系統在模型滾轉方向的慣量變化約為7‰,在偏航方向的慣量變化約為5‰,因此對振蕩頻率影響較小。

曲桿加工完成后,曲桿慣量基本不變。因此,為了更好地模擬真實尾旋運動的信息,機構設計時應嚴格控制曲桿慣量,如通過材料和工藝的選取以嚴格控制曲桿慣量所占模型的慣量比重。根據仿真結果,建議將曲桿和模型的滾轉慣量比值控制在10%以內。

4.3 摩擦力的影響

3自由度虛擬飛行試驗裝置中存在3個自由旋轉鉸,每個鉸聯接副之間存在固有的摩擦力矩。仿真分析了無摩擦、1倍摩擦和2倍摩擦力矩對尾旋試驗結果的影響,具體的對比結果見圖13。從圖中可以發現,不同摩擦力矩下模型在尾旋進入和穩定階段響應形態一致,摩擦力矩的影響并不明顯。

圖13 不同摩擦力矩仿真結果Fig.13 Simulation results with different friction moments

4.4 質心偏移的影響

對于動力學相似縮比模型,模型質心的Z方向位置一般不易調整;而且模型安裝時也可能會出現模型質心與萬向鉸旋轉中心有偏差的情況。因此通過設置模型質心在Z軸方向上與萬向鉸旋轉中心的偏移量來進行仿真,以分析質心偏移帶來的影響。

如圖14所示,模型重心與萬向鉸旋轉中心不重合時,在初始配平舵偏下模型不再平衡,略有低頭;給定舵面操縱信號后,仍能進入尾旋,但隨著偏移量增大,各狀態量的振蕩幅值相應增大,重心偏移0、3、5 mm情況下,迎角振蕩幅值分別為4°、10°、16°,側滑角振蕩幅值分別為5°、11°、18°;均值基本不變,振蕩周期略微增加。

圖14 不同質心偏移仿真結果Fig.14 Results of simulation with different offsets of center of mass

5 結 論

基于低速風洞3自由度虛擬飛行試驗裝置進行了某典型飛機尾旋的建模和仿真研究。

1) 該裝置可在水平風洞實現飛機尾旋的進入、發展和改出各階段的模擬,其中穩定尾旋的運動參數及模型姿態與立式風洞尾旋試驗結果較為吻合,可作為飛機尾旋風洞試驗研究的有力補充。

2) 曲桿慣量不影響模型尾旋運動的振蕩均值、旋轉速率等特征參數,但會延長模型進入尾旋過程,影響尾旋特征參數的振蕩幅值。建議將曲桿和模型的滾轉慣量比值控制在10%以內。

3) 摩擦力矩在一定范圍內對模型尾旋運動的主要特征參數沒有明顯影響。

4) 質心偏移會影響模型尾旋運動的振蕩幅值,對均值、旋轉速率影響不大。建議控制模型重心與萬向鉸旋轉中心的偏移量在3 mm以內。

猜你喜歡
飛機模型
一半模型
鷹醬想要“小飛機”
飛機失蹤
環球時報(2022-05-30)2022-05-30 15:16:57
國航引進第二架ARJ21飛機
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
“拼座飛機”迎風飛揚
當代陜西(2019年11期)2019-06-24 03:40:28
乘坐飛機
3D打印中的模型分割與打包
神奇飛機變變變
主站蜘蛛池模板: 日本午夜影院| 亚洲综合久久成人AV| 91视频99| 国产三级毛片| 青青热久免费精品视频6| 久久夜色精品| 精品综合久久久久久97| 亚洲婷婷六月| 欧洲精品视频在线观看| 55夜色66夜色国产精品视频| 国产AV毛片| 天天色天天操综合网| 国产精品美女自慰喷水| 一级黄色网站在线免费看| 久久国产精品麻豆系列| 亚洲91在线精品| 综合五月天网| 亚洲精品视频免费| 精品无码视频在线观看| 亚洲国产成人精品一二区| 亚洲精品国产综合99| 久久久久夜色精品波多野结衣| 欧美丝袜高跟鞋一区二区| 国产99在线| 亚洲区一区| 欧美一级一级做性视频| 国禁国产you女视频网站| 91青草视频| 国产精品私拍99pans大尺度| 中文字幕在线播放不卡| 亚洲欧美国产五月天综合| 午夜啪啪网| 国产sm重味一区二区三区| 婷婷亚洲视频| 久久久亚洲色| 国产福利观看| 欧美日韩精品综合在线一区| 97人妻精品专区久久久久| 亚洲浓毛av| 粉嫩国产白浆在线观看| 97视频在线观看免费视频| 欧美成人精品在线| 国产永久在线视频| 日韩亚洲高清一区二区| 熟妇丰满人妻av无码区| 亚洲综合二区| 国产精品尤物在线| 毛片视频网址| 国产成人精品男人的天堂下载| 亚洲色欲色欲www网| 999国产精品永久免费视频精品久久 | 一级毛片基地| 欧美成一级| 久久综合九色综合97婷婷| 伊人91在线| 亚洲美女视频一区| 丁香综合在线| 国产精品30p| 日韩在线1| 国产在线拍偷自揄拍精品| 亚洲六月丁香六月婷婷蜜芽| 99在线观看精品视频| 在线精品亚洲一区二区古装| 婷婷六月综合网| 色综合久久久久8天国| 亚洲男人天堂2020| 毛片免费高清免费| 99在线视频网站| 亚洲第一视频免费在线| 一本大道在线一本久道| 午夜视频在线观看免费网站| 成人午夜视频免费看欧美| 日韩视频免费| 91在线中文| 露脸真实国语乱在线观看| 伊人久综合| 欧美中文一区| 老司机午夜精品网站在线观看| 91成人免费观看在线观看| 91网站国产| 亚洲最新地址| 国产日本欧美亚洲精品视|