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

低速飛機(jī)分布式非線性氣動(dòng)模型建模方法

2012-11-06 09:08:49熊磊周洲
飛行力學(xué) 2012年5期
關(guān)鍵詞:飛機(jī)方法

熊磊, 周洲

(西北工業(yè)大學(xué) 無人機(jī)特種技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710065)

低速飛機(jī)分布式非線性氣動(dòng)模型建模方法

熊磊, 周洲

(西北工業(yè)大學(xué) 無人機(jī)特種技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710065)

針對(duì)低速類大展弦比飛機(jī)陣風(fēng)載荷分析與風(fēng)場(chǎng)飛行的仿真問題,提出了一種基于二維翼型氣動(dòng)特性的全機(jī)三維分布式非線性氣動(dòng)模型的建模方法。將采用該方法的計(jì)算結(jié)果與采用CFD的計(jì)算結(jié)果進(jìn)行了對(duì)比,證明了該方法可以精確模擬低速類大展弦比飛機(jī)的氣動(dòng)特性。以此為基礎(chǔ),進(jìn)行了無人機(jī)在1-Cosine型陣風(fēng)作用下的飛行仿真。結(jié)果表明,低速類大展弦比飛機(jī)在遇到陣風(fēng)時(shí),會(huì)出現(xiàn)與傳統(tǒng)問題完全不同的陣風(fēng)過載與飛行特性。

翼型; 分布式; 非線性氣動(dòng)模型; 陣風(fēng)載荷; 飛行仿真

引言

飛機(jī)在飛行時(shí)會(huì)受到來自大氣風(fēng)場(chǎng)的擾動(dòng)。用傳統(tǒng)方法[1-4]分析飛機(jī)遭遇陣風(fēng)擾動(dòng)時(shí),大都將其視為質(zhì)點(diǎn)。大型飛機(jī)由于受風(fēng)場(chǎng)影響更加顯著,因此采用了較為復(fù)雜的四點(diǎn)法或五點(diǎn)法。這些方法對(duì)傳統(tǒng)問題都有較好的適用性,但對(duì)于近年來出現(xiàn)的大展弦比、低翼載荷、低速度的新型飛機(jī),如太陽(yáng)能無人機(jī),則存在比較明顯的缺陷。此類飛機(jī)外形尺寸巨大,在風(fēng)場(chǎng)中會(huì)明顯受到風(fēng)速分布變化帶來的影響。此外,由于飛行速度低,在風(fēng)速足夠大的情況下,采用傳統(tǒng)方法將不再合適。例如陣風(fēng)風(fēng)向與飛機(jī)飛行方向夾角很大時(shí),飛機(jī)的相對(duì)迎角已遠(yuǎn)遠(yuǎn)超過其失速迎角,導(dǎo)致氣動(dòng)力發(fā)生非線性變化,且風(fēng)速會(huì)與飛行速度疊加,顯著增加當(dāng)?shù)貏?dòng)壓。以上這些問題都是傳統(tǒng)方法沒有或缺乏考慮的。

為了同時(shí)模擬分布式風(fēng)場(chǎng)和非線性氣動(dòng)力變化對(duì)飛機(jī)飛行特性帶來的影響,并在將來能夠加入對(duì)飛機(jī)彈性變形的考慮,本文結(jié)合文獻(xiàn)[5]中的思想,提出了一種基于二維翼型氣動(dòng)特性的全機(jī)三維分布式非線性氣動(dòng)模型建模方法。此方法的核心思想在于將飛機(jī)的氣動(dòng)特性進(jìn)行離散化分析,用飛機(jī)各氣動(dòng)部件的翼型氣動(dòng)特性對(duì)飛機(jī)的全機(jī)氣動(dòng)特性進(jìn)行合成。由于可以在很寬的迎角、速度及雷諾數(shù)范圍內(nèi)提供翼型的氣動(dòng)特性,且對(duì)于低速類飛機(jī)采用的大展弦比、小后掠角機(jī)翼又有著相對(duì)較高的精度,因此這種方法非常適合解決上述問題。

1 建模方法

1.1 飛機(jī)質(zhì)點(diǎn)運(yùn)動(dòng)向當(dāng)?shù)匾硇瓦\(yùn)動(dòng)的轉(zhuǎn)換

建模過程一共采用了3個(gè)坐標(biāo)系。重心處的機(jī)體坐標(biāo)系為模型的主坐標(biāo)系。當(dāng)?shù)貦C(jī)體坐標(biāo)系:該坐標(biāo)系遵從機(jī)體坐標(biāo)系原則,原點(diǎn)位于飛機(jī)任意位置;當(dāng)?shù)匾硇腕w坐標(biāo)系:該坐標(biāo)系兩個(gè)方向分別與當(dāng)?shù)匾硇偷亩S體坐標(biāo)系重合,而另外一個(gè)軸垂直于翼型平面指向展向。

在主坐標(biāo)系下,首先對(duì)飛機(jī)三維外形進(jìn)行離散。由于飛機(jī)外形為剛性,因此可將飛機(jī)的主要部件離散為有限個(gè)確定的坐標(biāo)點(diǎn)。將這些點(diǎn)的坐標(biāo)值與飛機(jī)重心處的速度和轉(zhuǎn)動(dòng)角速度結(jié)合,即可求出此點(diǎn)在當(dāng)?shù)貦C(jī)體坐標(biāo)系下的運(yùn)動(dòng)速度及轉(zhuǎn)動(dòng)角速度。設(shè)編號(hào)為k的坐標(biāo)點(diǎn)處的坐標(biāo)值為(xk,yk,zk),飛機(jī)重心在主坐標(biāo)系下的運(yùn)動(dòng)速度及轉(zhuǎn)動(dòng)角速度分別為vx,vy,vz,ωx,ωy,ωz,則此坐標(biāo)點(diǎn)處在當(dāng)?shù)伢w軸系下的運(yùn)動(dòng)速度及轉(zhuǎn)動(dòng)角速度分別為:

(1)

獲得當(dāng)?shù)伢w軸系下的速度與角速度后,再對(duì)當(dāng)?shù)匾硇偷陌惭b角、扭轉(zhuǎn)角和上反角進(jìn)行坐標(biāo)旋轉(zhuǎn),變換為當(dāng)?shù)匾硇腕w坐標(biāo)系。設(shè)此點(diǎn)處的機(jī)翼上反角為Γ,左機(jī)翼上反時(shí)取正;安裝角和扭轉(zhuǎn)角分別為i和ε,機(jī)翼前緣向上時(shí)取正,則變換矩陣A為:

(2)

1.2 二維翼型氣動(dòng)特性的表征

通過坐標(biāo)變換已經(jīng)將翼型的運(yùn)動(dòng)分離為翼型弦向及機(jī)翼展向兩個(gè)方向。對(duì)于大展弦比機(jī)翼可以忽略展向運(yùn)動(dòng),翼型的運(yùn)動(dòng)將只剩下平動(dòng)與轉(zhuǎn)動(dòng)兩部分。平動(dòng)部分可以通過兩個(gè)方向的平動(dòng)速度與當(dāng)?shù)貧饬?如陣風(fēng))疊加轉(zhuǎn)化成定常氣動(dòng)問題,而轉(zhuǎn)動(dòng)部分可將翼型的體軸坐標(biāo)原點(diǎn)設(shè)在50%弦長(zhǎng)處。由于前后弦長(zhǎng)相等,可近似認(rèn)為翼型繞此點(diǎn)的旋轉(zhuǎn)產(chǎn)生一對(duì)力偶,無作用力但有阻尼力矩。由于轉(zhuǎn)動(dòng)角速度一般非常小,因此可近似為線性值。

翼型平動(dòng)時(shí)的氣動(dòng)特性是本方法的核心。翼型在其二維體坐標(biāo)系下的氣動(dòng)特性主要靠弦向力、法向力與俯仰力矩3個(gè)因子來表征。其中,弦向力由于翼型外形呈狹長(zhǎng)型,且主要由摩擦及翼型前后壓差產(chǎn)生,因此造成其作用點(diǎn)變化范圍較小,可近似認(rèn)為作用于翼型弦線上。因而翼型力矩的變化可看作是法向力作用點(diǎn)變化所引起的。設(shè)翼型體軸系下法向力系數(shù)為CN,俯仰力矩系數(shù)為Cm,翼型當(dāng)?shù)叵议L(zhǎng)為cA,法向力作用點(diǎn)與力矩參考點(diǎn)的距離為xAF,則:

xAF=cACm/CN

(3)

使用翼型氣動(dòng)特性來反推全機(jī)氣動(dòng)特性的最大優(yōu)勢(shì)在于可采用分布式風(fēng)場(chǎng),并降低了飛行仿真中處理非線性氣動(dòng)導(dǎo)數(shù)的難度;可以對(duì)飛機(jī)各點(diǎn)處施加不同的風(fēng)速,這更接近真實(shí)物理過程,且理論上只要翼型數(shù)據(jù)準(zhǔn)確,飛機(jī)的各氣動(dòng)數(shù)據(jù)也應(yīng)是準(zhǔn)確的。本文中的翼型數(shù)據(jù)由CFD計(jì)算得到。理論上,采用風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)將更加合理,但在缺乏風(fēng)洞實(shí)驗(yàn)條件且只是定性分析問題時(shí),CFD的結(jié)果可以滿足需要。

翼型氣動(dòng)數(shù)據(jù)還可準(zhǔn)確模擬當(dāng)?shù)貋砹鞯乃俣茸兓?捕捉動(dòng)壓變化對(duì)氣動(dòng)特性的影響。這對(duì)在非均勻風(fēng)場(chǎng)下飛行的低速飛機(jī)來說有實(shí)際意義。而通過建立若干不同雷諾數(shù)對(duì)應(yīng)的翼型氣動(dòng)參數(shù),在其之間進(jìn)行插值則可以很好地模擬低速時(shí)的雷諾數(shù)效應(yīng)。

1.3 二維翼型氣動(dòng)特性向全機(jī)氣動(dòng)特性的轉(zhuǎn)換

在得到二維翼型的氣動(dòng)力與力作用點(diǎn)后,需要將其反變換回到全機(jī)的主坐標(biāo)系之下。氣動(dòng)力的變換矩陣為AT。力矩則需要考慮翼型法向力在翼型弦向方向的移動(dòng),其公式為:

(4)

對(duì)于大展弦比機(jī)翼,采用經(jīng)典的升力線理論進(jìn)行修正是合理的。文獻(xiàn)[6]提供了一種非常好的非線性氣動(dòng)力的升力線修正方法。飛機(jī)各主要部件之間的干擾修正需要針對(duì)對(duì)象的不同而采用合適的方法。本文采用文獻(xiàn)[7]提供的工程方法對(duì)翼面間的下洗和阻滯進(jìn)行修正。修正結(jié)果表明,對(duì)于低雷諾數(shù)飛機(jī)來說,修正值有些偏大。飛機(jī)非氣動(dòng)部件的氣動(dòng)力估算也按照文獻(xiàn)[7]中的對(duì)應(yīng)方法進(jìn)行修正。以上計(jì)算完成之后,對(duì)飛機(jī)各部件氣動(dòng)力進(jìn)行總成,即是飛機(jī)在給定飛行條件下的氣動(dòng)總力及力矩。在得到飛機(jī)質(zhì)心處的總力與力矩后,只需加入最基本形式的六自由度方程即可進(jìn)行飛行仿真。

2 算例分析

2.1 驗(yàn)證算例

此次模擬的對(duì)象為某型太陽(yáng)能無人機(jī),該無人機(jī)采用鴨式布局,翼載荷為5.54 kg/m2。計(jì)算高度為海平面,計(jì)算速度為10 m/s。本文方法的計(jì)算結(jié)果與CFD的計(jì)算結(jié)果對(duì)比如圖1~圖3所示。

圖1 升力系數(shù)隨迎角的變化曲線

圖2 阻力系數(shù)隨迎角的變化曲線

圖3 俯仰力矩系數(shù)隨迎角的變化曲線

從圖上可以看出,飛機(jī)的升力與力矩特性吻合非常好,而阻力值相對(duì)較小,這主要是前后翼干擾的工程修正法不夠精確所致,進(jìn)一步的工作應(yīng)采用風(fēng)洞或CFD數(shù)據(jù)進(jìn)行修正。橫航向氣動(dòng)特性的對(duì)比則證明由本文方法得到的側(cè)向力與滾轉(zhuǎn)力矩非常準(zhǔn)確,航向力矩則由于無法模擬大側(cè)滑時(shí)垂尾的氣動(dòng)效應(yīng)損失而偏保守,部分動(dòng)導(dǎo)數(shù)量級(jí)正確但存在差距,由于動(dòng)導(dǎo)數(shù)是由工程方法得到的,因此需要進(jìn)一步的分析。

2.2 1-Cosine型風(fēng)場(chǎng)飛行仿真及分析

在前文驗(yàn)證算例的基礎(chǔ)上,對(duì)某型太陽(yáng)能無人機(jī)進(jìn)行了1-Cosine型陣風(fēng)作用下的風(fēng)場(chǎng)飛行仿真。陣風(fēng)的速度型如下:

(5)

仿真共考慮了兩種典型的陣風(fēng)作用情況。第1種情況假設(shè)飛機(jī)整體進(jìn)入陣風(fēng),以模擬大尺度陣風(fēng)短時(shí)間作用于飛機(jī)上的過程。在此基礎(chǔ)上,考慮了兩種飛行方式。首先限制飛機(jī)的俯仰轉(zhuǎn)動(dòng),以觀察飛控系統(tǒng)控制律對(duì)過載的影響。圖4是過載的時(shí)間歷程曲線。一共模擬了5個(gè)陣風(fēng)作用角度。通過曲線對(duì)比可以看到,當(dāng)陣風(fēng)為垂直陣風(fēng)時(shí),過載增量并不大,因?yàn)榇藭r(shí)飛機(jī)已進(jìn)入深失速狀態(tài)。45°對(duì)應(yīng)的是飛機(jī)未進(jìn)入失速所受到的最大過載。60°時(shí)產(chǎn)生最大過載的原因在于飛機(jī)退出失速后與大速壓相互疊加。由于飛機(jī)翼載低,運(yùn)動(dòng)響應(yīng)顯著,因此失速前的最大過載處并非對(duì)應(yīng)出現(xiàn)在陣風(fēng)速度最大處,而是出現(xiàn)在陣風(fēng)速度增長(zhǎng)率最高處附近。而當(dāng)陣風(fēng)風(fēng)速增加速率減緩,導(dǎo)致飛機(jī)退出失速時(shí),陣風(fēng)速度卻達(dá)到最大,造成一個(gè)過載峰值。

作為對(duì)比,在相同參數(shù)下放開飛機(jī)的俯仰控制,此時(shí)飛機(jī)對(duì)陣風(fēng)自由響應(yīng),結(jié)果如圖5所示。由圖可以看到,在中小迎角時(shí),飛機(jī)的自由俯仰運(yùn)動(dòng)可以減小陣風(fēng)過載。這是因?yàn)轱w機(jī)在靜穩(wěn)定性作用下低頭,減小了來流迎角與陣風(fēng)過載。但當(dāng)陣風(fēng)夾角接近垂直時(shí),飛機(jī)同樣低頭減小來流角,但卻使飛機(jī)遇到了類似圖4中60°陣風(fēng)角時(shí)的情況,導(dǎo)致過載增加。因此在進(jìn)行陣風(fēng)過載抑制的控制系統(tǒng)設(shè)計(jì)時(shí),需要根據(jù)陣風(fēng)與飛機(jī)飛行速度方向綜合考慮。

圖4 過載曲線(有俯仰控制)

圖5 過載曲線(無俯仰控制)

第2種情況假設(shè)飛機(jī)沿飛行方向先后進(jìn)入一個(gè)垂直陣風(fēng)。此情況主要研究風(fēng)場(chǎng)隨時(shí)間的變化(即非凍結(jié)場(chǎng))及飛機(jī)縱向尺寸對(duì)飛機(jī)過載的影響。為了模擬風(fēng)場(chǎng)的時(shí)間變化,加入時(shí)間延遲tD=d/vD,式中,d為飛機(jī)上某點(diǎn)最初進(jìn)入陣風(fēng)點(diǎn)的距離,vD為飛機(jī)在該方向的飛行速度與陣風(fēng)速度的疊加值。為了研究飛機(jī)的基本特性對(duì)其風(fēng)場(chǎng)飛行特性的影響,因此賦予其完全的自由度。圖6所示的是相同垂直陣風(fēng)作用下沒有考慮風(fēng)場(chǎng)時(shí)變與考慮了風(fēng)場(chǎng)時(shí)變的過載響應(yīng)曲線對(duì)比。

圖6 過載對(duì)比曲線(無俯仰控制)

由圖可以看到,飛機(jī)過載會(huì)有一個(gè)明顯的延遲。由于尺寸和飛行速度低,導(dǎo)致飛機(jī)鴨翼首先進(jìn)入陣風(fēng),并使飛機(jī)在初期上仰。兩次仿真的最大過載點(diǎn)特征相似,這是飛機(jī)機(jī)翼產(chǎn)生最大升力的標(biāo)志。但加入風(fēng)場(chǎng)時(shí)變后的過載值較小,這是鴨翼與機(jī)翼的最大升力點(diǎn)相互錯(cuò)開所造成的。由于飛機(jī)機(jī)翼進(jìn)入最大陣風(fēng)點(diǎn)時(shí)鴨翼處的風(fēng)速已經(jīng)非常小,因此低頭力矩更大并造成飛機(jī)的低頭運(yùn)動(dòng)更加劇烈,飛機(jī)在陣風(fēng)結(jié)束時(shí)俯沖角較非時(shí)變風(fēng)場(chǎng)大很多。大展弦比飛機(jī)的縱向轉(zhuǎn)動(dòng)慣量一般會(huì)比其他兩個(gè)轉(zhuǎn)動(dòng)慣量小一個(gè)量級(jí),這在很大程度上惡化了縱向響應(yīng)問題。

綜上所述,風(fēng)場(chǎng)隨時(shí)間的變化對(duì)這類飛機(jī)的運(yùn)動(dòng)響應(yīng)有著顯著的影響,因此在研究此類問題特別是紊流場(chǎng)這種典型分布式風(fēng)場(chǎng)時(shí),必須考慮風(fēng)場(chǎng)隨時(shí)間變化對(duì)飛機(jī)飛行造成的影響。

3 結(jié)束語

本文采用分布式非線性氣動(dòng)模型對(duì)低速類大展弦比飛機(jī)的氣動(dòng)特性進(jìn)行了模擬。結(jié)果表明,該方法原理正確,精度可初步滿足工程設(shè)計(jì)和飛行仿真的需求,并具備較大的提升空間。由于該方法在非線性氣動(dòng)特性模擬、計(jì)算效率上具有較大優(yōu)勢(shì),并具備耦合結(jié)構(gòu)模型、模擬飛機(jī)彈性變形的潛力,因此非常適合應(yīng)用于低速類大展弦比飛機(jī)的飛行仿真領(lǐng)域。對(duì)太陽(yáng)能無人機(jī)1-Cosine型陣風(fēng)響應(yīng)的仿真結(jié)果表明,低速大展弦比飛機(jī)的過載產(chǎn)生機(jī)理、產(chǎn)生過程、大小均與常規(guī)飛機(jī)不同,超過傳統(tǒng)方法的適用范圍,并嚴(yán)重影響到飛機(jī)的飛行安全。此外,風(fēng)場(chǎng)的時(shí)變特性與飛機(jī)本身的尺寸對(duì)飛機(jī)飛行特性的影響非常顯著,在飛行仿真中必須考慮分布式時(shí)變風(fēng)場(chǎng)的影響。

[1] 飛機(jī)設(shè)計(jì)手冊(cè)總編委會(huì).飛機(jī)設(shè)計(jì)手冊(cè)[M].北京:航空工業(yè)出版社,1994.

[2] 國(guó)防科學(xué)技術(shù)工業(yè)委員會(huì).GJB 67.2-85 軍用飛機(jī)強(qiáng)度和剛度規(guī)范——飛行載荷[S].北京:國(guó)防科學(xué)技術(shù)工業(yè)委員會(huì),1985.

[3] Frederic M Hoblit.Gust loads on aircraft:concepts and applications[M].Washington,D C:American Institute of Aeronautics and Astronautics Inc,1988.

[4] 肖業(yè)論,金長(zhǎng)江.大氣擾動(dòng)中的飛行原理[M].北京:國(guó)防工業(yè)出版社,1993.

[5] 何植岱,高浩.高等飛行動(dòng)力學(xué)[M].西安:西北工業(yè)大學(xué)出版社,1990.

[6] 吳子牛.空氣動(dòng)力學(xué)[M]. 北京:清華大學(xué)出版社,2007.

[7] 嚴(yán)恒元,陳勁松,馮亞南,等.飛行器氣動(dòng)特性分析與工程計(jì)算[M].西安:西北工業(yè)大學(xué)出版社,1990.

(編輯:姚妙慧)

Amethodfordistributingnonlinearaerodynamicsmodeloflow-speedaircraft

XIONG Lei, ZHOU Zhou

(National Key Laboratory of Science and Technology on UAV, NWPU, Xi’an 710065, China)

A method for the whole aircraft 3-dimemsional distributed nonlinear aerodynamics model based on the aerodynamics of 2-dimensional airfoils has been established in order to calculate the wind load and simulate the flight process of low-speed aircraft with high-aspect-ratio in wind field. This method has already been proved reasonable by modeling a solar aircraft and comparing to CFD results. On the basis of this, a flight simulation through a 1-Cosine gust wind has also been done, from which it was found that the low speed aircraft with high aspect-ratio will have completely different gust load and flight characteristics when encountering gust wind.

airfoil; distributing; nonlinear aerodynamics model; wind load; flight simulation

V211.4

A

1002-0853(2012)05-0398-04

2011-11-23;

2012-04-16

熊磊(1985-),男,江西南昌人,碩士研究生,主要研究方向?yàn)轱w行器設(shè)計(jì)、飛行力學(xué)、計(jì)算流體力學(xué);

周洲(1966-),女,湖南長(zhǎng)沙人,教授,博士生導(dǎo)師,主要研究方向?yàn)闊o人機(jī)總體、氣動(dòng)布局設(shè)計(jì)。

猜你喜歡
飛機(jī)方法
鷹醬想要“小飛機(jī)”
飛機(jī)失蹤
國(guó)航引進(jìn)第二架ARJ21飛機(jī)
“拼座飛機(jī)”迎風(fēng)飛揚(yáng)
學(xué)習(xí)方法
乘坐飛機(jī)
神奇飛機(jī)變變變
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
主站蜘蛛池模板: 国产综合另类小说色区色噜噜| 国产毛片不卡| 亚洲无线观看| 国产一级片网址| 国产精品人人做人人爽人人添| 91久久偷偷做嫩草影院| 国产69精品久久| 国产午夜小视频| 色综合国产| 日韩东京热无码人妻| 99热这里只有精品2| 四虎永久在线| 婷婷色丁香综合激情| 人妻中文字幕无码久久一区| 伊人激情综合网| 99久久亚洲综合精品TS| 3344在线观看无码| 亚洲精品色AV无码看| 日本高清免费不卡视频| 婷婷综合缴情亚洲五月伊| 99视频精品全国免费品| 亚洲精品第一页不卡| 亚洲AV一二三区无码AV蜜桃| 丁香亚洲综合五月天婷婷| 色视频国产| 国产成人综合在线观看| 久久九九热视频| 国产小视频免费| 国产真实二区一区在线亚洲| 幺女国产一级毛片| 亚洲精品卡2卡3卡4卡5卡区| 91小视频版在线观看www| 老司国产精品视频| 制服丝袜 91视频| 日本欧美午夜| 高清欧美性猛交XXXX黑人猛交| 国产女人18水真多毛片18精品| 欧美成人午夜视频| 国产成人精品一区二区三在线观看| 久久男人视频| 91人妻在线视频| 亚洲人视频在线观看| 亚洲色欲色欲www在线观看| 日韩成人在线视频| 国产伦精品一区二区三区视频优播| 国产精品对白刺激| 国产超薄肉色丝袜网站| 国产精品一区二区无码免费看片| 国产精品成人久久| 欧美日韩中文字幕在线| 韩国自拍偷自拍亚洲精品| 亚洲最黄视频| 看你懂的巨臀中文字幕一区二区| 久久久噜噜噜久久中文字幕色伊伊 | 成人在线不卡| 免费国产不卡午夜福在线观看| 国产91全国探花系列在线播放| 视频二区亚洲精品| 亚洲国产黄色| 狠狠色综合网| 国产第二十一页| 亚洲第一福利视频导航| 久久精品丝袜高跟鞋| 亚洲动漫h| 女人18毛片一级毛片在线 | 无码人妻免费| 亚洲精品成人7777在线观看| 成年午夜精品久久精品| 欧美日韩国产在线人| 国产精品自拍露脸视频| 青青草一区| 免费va国产在线观看| 亚洲精品国产首次亮相| 日韩人妻精品一区| 国产理论精品| 伊人久久精品无码麻豆精品| 欧美天堂在线| 国产打屁股免费区网站| 香蕉在线视频网站| 国产在线观看99| 亚洲无码视频一区二区三区 | 色成人亚洲|