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

高超聲速飛行器剛體/彈性體耦合動(dòng)力學(xué)建模

2012-06-22 06:59:56李惠峰
關(guān)鍵詞:模態(tài)振動(dòng)結(jié)構(gòu)

李惠峰 肖 進(jìn) 張 冉

(北京航空航天大學(xué) 宇航學(xué)院,北京 100191)

高超聲速飛行器剛體/彈性體耦合動(dòng)力學(xué)建模

李惠峰 肖 進(jìn) 張 冉

(北京航空航天大學(xué) 宇航學(xué)院,北京 100191)

高超聲速飛行器廣泛采用升力體、乘波體等氣動(dòng)布局和輕質(zhì)復(fù)合材料、薄壁結(jié)構(gòu)等,導(dǎo)致結(jié)構(gòu)振動(dòng)與剛體運(yùn)動(dòng)頻率非常接近,給飛行器制導(dǎo)控制系統(tǒng)設(shè)計(jì)帶來(lái)了巨大挑戰(zhàn).針對(duì)該類(lèi)飛行器的特點(diǎn),考慮結(jié)構(gòu)的橫向位移,將機(jī)身前后體簡(jiǎn)化為于質(zhì)心處固聯(lián)的2根懸臂梁,并從統(tǒng)一的能量觀點(diǎn)出發(fā),基于拉格朗日方程與虛功原理,在縱向平面推導(dǎo)出適合高超聲速飛行器的剛體/彈性體耦合動(dòng)力學(xué)模型.通過(guò)對(duì)比耦合模型與傳統(tǒng)剛體模型的極點(diǎn)分布情況,發(fā)現(xiàn)結(jié)構(gòu)振動(dòng)與剛體短周期模態(tài)緊密耦合,離心力的引入影響了高度與長(zhǎng)周期模態(tài),對(duì)高超聲速飛行器航跡運(yùn)動(dòng)的作用不可忽視.最后分析了飛行速度與結(jié)構(gòu)阻尼變化對(duì)耦合模型動(dòng)態(tài)性能的影響.結(jié)果證明飛行速度對(duì)剛體運(yùn)動(dòng)模態(tài)影響顯著,而結(jié)構(gòu)阻尼的變化主要改變彈性模態(tài).

高超聲速;耦合模型;結(jié)構(gòu)振動(dòng);拉格朗日方程

隨著高超聲速飛行器的出現(xiàn),為節(jié)省燃料廣泛采用輕質(zhì)的復(fù)合材料,設(shè)計(jì)細(xì)長(zhǎng)體、升力體等氣動(dòng)布局,這些因素導(dǎo)致結(jié)構(gòu)固有振動(dòng)頻率的降低,于是必須考慮剛體和彈性體的慣性耦合現(xiàn)象[1].特別是以吸氣式超燃沖壓發(fā)動(dòng)機(jī)為動(dòng)力的X-43和X-51,發(fā)動(dòng)機(jī)對(duì)飛行條件以及氣動(dòng)外形提出的要求非常苛刻,在高馬赫數(shù)下機(jī)體的彈性振動(dòng)直接影響到激波角,進(jìn)而改變了發(fā)動(dòng)機(jī)的進(jìn)氣情況.

以協(xié)和號(hào)超音速客機(jī)為例首次提出剛體與彈性振動(dòng)的耦合建模[2],但是耦合僅體現(xiàn)在定常與非定常氣動(dòng)力的相互影響.在此基礎(chǔ)上,文獻(xiàn)[3]推導(dǎo)出剛體/彈性體的動(dòng)力學(xué)耦合一般模型,稱(chēng)為強(qiáng)耦合真實(shí)模型(TMH,Truth Model with Heave coupling).但TMH中未考慮高馬赫數(shù)下地球離心力對(duì)飛行器運(yùn)動(dòng)的影響,且未將推力加入到機(jī)體軸向的非保守力系,不能真實(shí)地體現(xiàn)高超聲速飛行器的特點(diǎn).文獻(xiàn)[4]將TMH中相對(duì)較弱的耦合項(xiàng)去掉,得到基于控制的簡(jiǎn)化模型(TM,Truth Model).文獻(xiàn)[5]將尖頭體高超聲速飛行器機(jī)體前緣的結(jié)構(gòu)振動(dòng)速度與攻角狀態(tài)線性疊加,主要是為了驗(yàn)證結(jié)構(gòu)振動(dòng)對(duì)發(fā)動(dòng)機(jī)進(jìn)氣的影響,但其忽略了對(duì)其他飛行狀態(tài)量的影響以及剛體運(yùn)動(dòng)對(duì)結(jié)構(gòu)振動(dòng)的反作用.為了更真實(shí)地體現(xiàn)飛行器結(jié)構(gòu)振動(dòng)模態(tài),文獻(xiàn)[6]將機(jī)身簡(jiǎn)化為一根自由-自由梁,取前三階固有振型,結(jié)果顯示慣性耦合項(xiàng)系數(shù)為零,無(wú)法體現(xiàn)出結(jié)構(gòu)與剛體運(yùn)動(dòng)的強(qiáng)耦合現(xiàn)象.綜上所述,國(guó)外已進(jìn)行過(guò)耦合建模的相關(guān)研究,但研究對(duì)象及內(nèi)容各有不同.由于我國(guó)在高超聲速飛行器的理論研究和實(shí)際應(yīng)用方面相對(duì)落后,本文研究目的是在充分借鑒國(guó)外公開(kāi)發(fā)表的文獻(xiàn)資料基礎(chǔ)上,結(jié)合吸氣式高超聲速飛行器的特點(diǎn),推導(dǎo)出較真實(shí)的動(dòng)力學(xué)耦合模型,以便為未來(lái)針對(duì)該類(lèi)型飛行器的先進(jìn)制導(dǎo)和控制系統(tǒng)理論研究和仿真驗(yàn)證以及進(jìn)行結(jié)構(gòu)/制導(dǎo)控制一體化設(shè)計(jì)打下良好基礎(chǔ).

本文首先建立機(jī)體的結(jié)構(gòu)振動(dòng)模型,將機(jī)體簡(jiǎn)化為2根于質(zhì)心處固連的懸臂梁,分別求出前四階的固有振型和自然振動(dòng)頻率;然后建立剛體/彈性體的耦合矢量模型,結(jié)合吸氣式高超聲速飛行器的特點(diǎn),基于拉格朗日方程推導(dǎo)出耦合動(dòng)力學(xué)模型;最后通過(guò)將耦合非線性模型線性化,求取平衡點(diǎn),并在此基礎(chǔ)上分析極點(diǎn)分布情況,總結(jié)了結(jié)構(gòu)振動(dòng)與剛體運(yùn)動(dòng)的耦合原因以及離心力、飛行速度與結(jié)構(gòu)阻尼對(duì)耦合模型動(dòng)態(tài)性能的影響.

1 機(jī)身振動(dòng)簡(jiǎn)化模型

在飛行器設(shè)計(jì)的初級(jí)分析階段,為了減少仿真計(jì)算量的同時(shí)又不失真實(shí)性,通常將全機(jī)振動(dòng)模型簡(jiǎn)化為彈性梁.本文以類(lèi)X-43A飛行器為研究對(duì)象,由于此飛行器長(zhǎng)度遠(yuǎn)大于厚度,機(jī)身近似細(xì)長(zhǎng)體,并考慮到推進(jìn)系統(tǒng)對(duì)前體激波的形成與后體外膨脹段的外形要求,機(jī)體設(shè)計(jì)成中間粗、兩端細(xì),容積和大部分質(zhì)量集中于質(zhì)心周?chē)?因此,將機(jī)身前體和后體分別簡(jiǎn)化為懸臂梁,于質(zhì)心處固連,是較合理的假設(shè),并且能體現(xiàn)出剛體與彈性體之間的強(qiáng)耦合[5].

結(jié)構(gòu)的變形分為橫向、軸向位移與扭轉(zhuǎn)變形,為簡(jiǎn)化起見(jiàn),本文只研究變形最明顯的橫向位移.

基于以上假設(shè),振動(dòng)位移w(x,t)與機(jī)體軸向坐標(biāo)和時(shí)間有關(guān),由梁的彎曲理論,無(wú)阻尼自然振動(dòng)方程如下:

以前體懸臂梁為例,通過(guò)分離變量法,并結(jié)合邊界條件可求出結(jié)構(gòu)固有振型φ與自然振動(dòng)頻率ω的解析解:

其中,i為振型的階數(shù).

圖1是機(jī)體前懸臂梁的前四階振型,此振型值體現(xiàn)在下文中建立的剛體/彈性體慣性耦合項(xiàng)系數(shù)中.

對(duì)應(yīng)的自然振動(dòng)頻率分別為

通常,固有頻率越低,與剛體運(yùn)動(dòng)頻率越接近,故耦合越嚴(yán)重,本文對(duì)于前后體懸臂梁均取第1階振型與固有振動(dòng)頻率.

圖1 機(jī)身前端簡(jiǎn)化懸臂梁前四階振型

2 剛體/彈性體矢量模型建立

剛體/彈性體耦合建模的重點(diǎn)在于矢量模型的建立.不同于以往剛體假設(shè)的是,真實(shí)彈性體飛行器上各點(diǎn)與質(zhì)心的距離不是一成不變的,而是在其附近存在著微幅振動(dòng).本文假設(shè)只存在結(jié)構(gòu)的橫向位移,而忽略影響較小的軸向位移.矢量模型如圖2所示.

假設(shè)rO為從慣性系原點(diǎn)On指向機(jī)體系原點(diǎn)(質(zhì)心)O的矢量,如圖2,向量之間存在關(guān)系:

其中,O為飛行器質(zhì)心;P為剛體飛行器縱向平面上任一點(diǎn).現(xiàn)假設(shè)該飛行器伴隨有結(jié)構(gòu)彈性振動(dòng)位移:

其中,ηz,i為第 i階結(jié)構(gòu)橫向位移變量.

ρ=xxb+zzb表示剛體模型上P點(diǎn)相對(duì)于O點(diǎn)的位移矢量,所以有rP/O=d+ρ.于是,式(4)可變成 rP=rO+d+ρ.

圖2 結(jié)構(gòu)振動(dòng)/剛體運(yùn)動(dòng)耦合矢量模型

P點(diǎn)的速度矢量:

在縱向平面內(nèi),慣性系與機(jī)體系存在如下轉(zhuǎn)換關(guān)系:

其中,θ為俯仰角.

將慣性系下的質(zhì)心矢量投影到機(jī)體系下得

于是,有

得到機(jī)體系下質(zhì)心的速度矢量為

其中,U,W分別為飛行速度沿機(jī)體軸 xb,zb的分量.

3 剛體/彈性體耦合動(dòng)力學(xué)模型

結(jié)合以上的矢量定義,將整個(gè)飛行器視為多質(zhì)點(diǎn)系,只考慮縱向平面的位移,并利用振型的正交性可求出總動(dòng)能如下:

其中,Iyy為慣性矩;λz,i,ψxz,i為慣性耦合項(xiàng)系數(shù).

總勢(shì)能為重力勢(shì)能與彈性勢(shì)能之和,即

其中,m 為飛行器質(zhì)量;ωz,i為結(jié)構(gòu)第 i階自然振動(dòng)頻率.

根據(jù)虛功原理:

式中,下標(biāo)f,a分別代表飛行器前體與后體;T為發(fā)動(dòng)機(jī)推力;-2ζωη·為結(jié)構(gòu)振動(dòng)的阻尼力,ζ為結(jié)構(gòu)振動(dòng)阻尼;N為引起結(jié)構(gòu)振動(dòng)的載荷力;Fx,F(xiàn)z分別為機(jī)體系下沿機(jī)體軸向和法向的合力;MA為俯仰力矩.通過(guò)機(jī)體系與氣流系的坐標(biāo)轉(zhuǎn)換,可得到升力(L)和阻力(D)的表達(dá)式:

結(jié)合拉格朗日函數(shù)和虛功表達(dá)式進(jìn)行拉格朗日方程求解,可推導(dǎo)出高超聲速飛行器剛體/彈性體耦合的動(dòng)力學(xué)方程組:

式中,慣性耦合項(xiàng)系數(shù)求解如下:

以上是機(jī)體系下的動(dòng)力學(xué)方程,通常所需的狀態(tài)量攻角(α)和速度(V)由下式表示:

通常情況下,低速飛行器的重力加速度求解只與高度有關(guān),但對(duì)于高超聲速飛行器,還必須考慮離心力對(duì)其的影響,所以經(jīng)簡(jiǎn)化的重力加速度求解公式為

其中,G為萬(wàn)有引力常數(shù);M為地球質(zhì)量;r為從地心到飛行器質(zhì)心的距離.

由于本文只在縱向平面建模,即假定橫側(cè)向狀態(tài)量均為0,使得與之相關(guān)的哥氏力在縱向平面的投影亦為0.并且,在高馬赫數(shù)下哥氏力相對(duì)于離心力為小量,因此,建模中只考慮離心力對(duì)耦合模型的作用.

4 模型分析

基于以上耦合模型,參考文獻(xiàn)[4]中的氣動(dòng)力模型及相關(guān)結(jié)構(gòu)、外形參數(shù)如下:

使用MATLAB中trim函數(shù)在固定高度與速度下求取平衡工作點(diǎn),見(jiàn)表1.其中2個(gè)控制輸入:發(fā)動(dòng)機(jī)節(jié)流閥值β和俯仰舵偏角δe.

由表1數(shù)據(jù)可看出,攻角越大,發(fā)動(dòng)機(jī)節(jié)流閥值越小,即省燃料;而結(jié)構(gòu)振動(dòng)幅值卻與之成正比,當(dāng)進(jìn)行攻角優(yōu)化時(shí)需折中考慮.

4.1 耦合模態(tài)分析

本文在高度為25 908 m,馬赫數(shù)為8平衡點(diǎn)下將耦合模型線性化并求解其特征根,同時(shí),為了體現(xiàn)出其特點(diǎn),求出了相同條件下純剛體運(yùn)動(dòng)模型的極點(diǎn)分布,以便對(duì)比.表2為耦合模型與剛體模型的極點(diǎn)分布情況.

表1 不同高度與速度下的平衡點(diǎn)

表2 極點(diǎn)分布對(duì)比

通過(guò)分析表2的數(shù)據(jù),可得出以下結(jié)論:

1)與經(jīng)典剛體模型不同,耦合模型由于加入了2個(gè)彈性振動(dòng)方程,總共出現(xiàn)9個(gè)特征根,可分為高度模態(tài)、長(zhǎng)周期模態(tài)、短周期模態(tài)和一對(duì)穩(wěn)定的氣動(dòng)彈性模態(tài).各運(yùn)動(dòng)參數(shù)隨時(shí)間的變化是上述4種模態(tài)的疊加.

2)與剛體模型的極點(diǎn)分布情況相比,二者的高度模態(tài)和長(zhǎng)周期模態(tài)非常接近,說(shuō)明巡航飛行中結(jié)構(gòu)振動(dòng)對(duì)高度、速度和航跡傾角的影響很小.但由于引入了結(jié)構(gòu)振動(dòng),剛體運(yùn)動(dòng)的2個(gè)短周期實(shí)根變成了2對(duì)共軛復(fù)根,說(shuō)明由于結(jié)構(gòu)的振動(dòng)導(dǎo)致了攻角和俯仰角速度的震蕩現(xiàn)象,同時(shí)結(jié)構(gòu)的振動(dòng)模態(tài)與短周期模態(tài)產(chǎn)生了耦合.由此可知,該耦合模型主要改變了剛體運(yùn)動(dòng)的短周期模態(tài),同時(shí),攻角等短周期狀態(tài)量的改變也嚴(yán)重影響了結(jié)構(gòu)的振動(dòng)特性.

由圖3所示,在相同平衡點(diǎn)處將耦合模型線性化,求出極點(diǎn)分布,并與完全去除慣性耦合項(xiàng)的模型極點(diǎn)分布對(duì)比.容易看出,除了代表高度模態(tài)和長(zhǎng)周期模態(tài)的中間3點(diǎn)基本重合,其余極點(diǎn)位置相距較遠(yuǎn),形象地表示出剛體/彈性體耦合運(yùn)動(dòng)模型中各模態(tài)的強(qiáng)耦合現(xiàn)象.

圖3 耦合/解耦模型極點(diǎn)分布對(duì)比

4.2 離心力對(duì)模態(tài)的影響

表3分別給出了耦合模型中有無(wú)離心力的極點(diǎn)分布情況.

表3 離心力對(duì)極點(diǎn)分布的影響

從表3中數(shù)據(jù)能看出,離心力對(duì)高度模態(tài)的影響最大,達(dá)到97.8%,直觀上說(shuō)是減緩了高度狀態(tài)量的發(fā)散速度,其次是長(zhǎng)周期模態(tài)3.3%,而對(duì)短周期模態(tài)和氣動(dòng)彈性模態(tài)的影響可以忽略.由此得知,雖然離心力對(duì)姿態(tài)狀態(tài)量與結(jié)構(gòu)振動(dòng)的影響非常小,在控制器設(shè)計(jì)時(shí)可以不予考慮,但是對(duì)于高超聲速剛體/彈性體耦合模型的制導(dǎo)控制一體化設(shè)計(jì)時(shí)不能簡(jiǎn)單省略.

4.3 參數(shù)變化對(duì)模態(tài)的影響

由圖4可見(jiàn),在同一高度下,飛行馬赫數(shù)從8變化到10的過(guò)程中,短周期的2對(duì)共軛復(fù)根逐漸向虛軸靠攏,說(shuō)明姿態(tài)角狀態(tài)量的衰減速率隨著馬赫數(shù)的增大而減小,且趨于臨界穩(wěn)定;當(dāng)速度繼續(xù)增大,其振蕩頻率隨之升高.從以上推導(dǎo)的耦合動(dòng)力學(xué)方程中也能看出,馬赫數(shù)的變化對(duì)飛行姿態(tài)和彈性振動(dòng)的影響非常明顯.

圖4 速度變化對(duì)模態(tài)的影響

圖5 結(jié)構(gòu)阻尼對(duì)模態(tài)的影響

如圖5所示,當(dāng)結(jié)構(gòu)阻尼從0.1~0.7變化時(shí)氣動(dòng)彈性模態(tài)的衰減速度加快、振動(dòng)頻率明顯降低,而對(duì)其它模態(tài)的影響相對(duì)較小.

5 結(jié)論

本文基于拉格朗日方程與虛功原理推導(dǎo)了高超聲速飛行器的剛體/彈性體耦合動(dòng)力學(xué)模型,分析得到如下結(jié)論:

1)結(jié)構(gòu)振動(dòng)與剛體短周期模態(tài)的相互影響非常明顯,對(duì)該類(lèi)型飛行器控制系統(tǒng)設(shè)計(jì)時(shí)需要著重考慮;

2)離心力的引入主要改變了高度和長(zhǎng)周期模態(tài),使得高度模態(tài)的發(fā)散速度減慢;

3)隨著飛行速度的增加,各極點(diǎn)值向零靠攏,整個(gè)系統(tǒng)趨于臨界穩(wěn)定,結(jié)構(gòu)阻尼的變化主要影響彈性模態(tài).該模型能真實(shí)地反映高速飛行器的運(yùn)動(dòng)特性,體現(xiàn)出了強(qiáng)耦合、強(qiáng)非線性的特點(diǎn).

需要說(shuō)明的是,本文主要考慮減少計(jì)算量的同時(shí)方便分析該模型,因此,結(jié)構(gòu)載荷力中未引入因結(jié)構(gòu)振動(dòng)產(chǎn)生的非定常氣動(dòng)力,但在實(shí)際工程應(yīng)用中高馬赫數(shù)下的非定常氣動(dòng)力不可忽視,這部分工作可在以后開(kāi)展.

(References)

[1]楊超,許赟,謝長(zhǎng)川.高超聲速飛行器氣動(dòng)彈性力學(xué)研究綜述[J].航空學(xué)報(bào),2010,31(1):1-11

Yang Chao,Xu Yun,Xie Changchuan.Review of studies on aeroelasticity of hypersonic vehicles[J].ACTA Aeronautica et Astronautica sinica,2010,31(1):1-11(in Chinese)

[2]Waszak M R,Schmit D K.On the flight dynamics of aeroelastic vehicles[R].AIAA-1986-2077-388,1986

[3]Bolender M A,Doman D B.A non-linear model for the longitudinal dynamics of a hypersonic air-breathing vehicle[R].AIAA-2005-6255,2005

[4]Parker J T,Bolender M A,Doman D B.Control-oriented modeling of an air-breathing hypersonic vehicle[R].AIAA-27830-611,2007

[5]Clark A D,Mirmirani M D,et al.An aero-propulsion integrated elastic model of a generic air-breathing hypersonic vehicle[R].AIAA-2006-6560,2006

[6]Rodriguez A A,Dickeson J J,Sridharan S,et al.Control relevant modeling,analysis,and design for scramjet-powered hypersonic vehicles[R].AIAA-2009-7287-980,2009

Hypersonic vehicle rigid/elastic coupled dynamic modeling

Li Huifeng Xiao Jin Zhang Ran
(School of Astronautics,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)

Extensive use of lifting body,wave-rider aerodynamic layout and other composite materials,thin-walled structures,lead the frequency of structural vibration and rigid body motion very close.This yields a great challenge to the vehicle control system design.Taking into account the transverse displacement,a specific simplification which treats the vehicle body as two mass center fixed cantilever beams was adopted.And the hypersonic aircraft's rigid body/elastic coupling model was derived based on the principle of virtual work and Lagrange equations.After comparing the open-looped poles of new coupling model and traditional rigid body model,a conclusion that the short-period vibration mode and the structures vibration mode tightly coupled with each other was made.The centrifugal force affected the height and long-period mode,and the effect on the flight path can not be ignored.Finally,changes in flight speed and structural damping on the dynamic performance of coupled model were analyzed.Results show that flight speeds significantly affect the rigid body motion modes,while the structural damping mainly changes in elastic mode.

hypersonic;coupling model;structure vibration;Lagrange equation

V 212.1

A

1001-5965(2012)02-0160-06

2010-10-14;< class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間:

時(shí)間:2012-02-21 11:46;

CNKI:11-2625/V.20120221.1146.010

www.cnki.net/kcms/detail/11.2625.V.20120221.1146.010.html

李惠峰(1970-),女,陜西蒲城人,副教授,lihuifeng@buaa.edu.cn.

(編 輯:趙海容)

猜你喜歡
模態(tài)振動(dòng)結(jié)構(gòu)
振動(dòng)的思考
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
振動(dòng)與頻率
論結(jié)構(gòu)
中立型Emden-Fowler微分方程的振動(dòng)性
論《日出》的結(jié)構(gòu)
國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長(zhǎng)
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
UF6振動(dòng)激發(fā)態(tài)分子的振動(dòng)-振動(dòng)馳豫
主站蜘蛛池模板: 亚洲天堂视频网| 538国产视频| v天堂中文在线| 呦系列视频一区二区三区| 国产精品视频免费网站| 亚洲精品视频网| 精品伊人久久大香线蕉网站| 亚洲成人一区二区三区| 国产精品尹人在线观看| 四虎AV麻豆| 久久精品日日躁夜夜躁欧美| 国国产a国产片免费麻豆| 91丝袜乱伦| 精品视频免费在线| 国产美女在线观看| 国产高清在线丝袜精品一区| 国产95在线 | 国产日韩欧美中文| 亚洲天堂久久久| 午夜国产精品视频| 丁香婷婷激情网| 老司机久久精品视频| 亚洲国产日韩视频观看| 国产91成人| 波多野结衣无码AV在线| 毛片免费试看| 少妇高潮惨叫久久久久久| 中文字幕欧美日韩| 97久久精品人人做人人爽| 欧美成人看片一区二区三区 | 欧美激情视频二区| lhav亚洲精品| 国产永久在线观看| 婷婷99视频精品全部在线观看| 国产黄色片在线看| jizz国产视频| 欧美乱妇高清无乱码免费| 国产成人高精品免费视频| 日韩专区第一页| 伊人激情综合| 狼友av永久网站免费观看| 日本一区高清| Jizz国产色系免费| 国产女人喷水视频| 国产香蕉国产精品偷在线观看| 国产精品毛片一区| 亚洲第一视频网| 国产另类乱子伦精品免费女| 国产欧美成人不卡视频| 国产精品思思热在线| 成人综合久久综合| 国内嫩模私拍精品视频| 就去色综合| 国产一区二区色淫影院| 激情视频综合网| 在线观看免费黄色网址| 重口调教一区二区视频| 91口爆吞精国产对白第三集| 亚洲成人在线免费| 欧美日韩午夜| 国产精品天干天干在线观看| 亚洲精品成人7777在线观看| 国产微拍一区二区三区四区| 亚洲视频免费播放| 日韩精品无码一级毛片免费| 不卡视频国产| 曰韩免费无码AV一区二区| 欧美激情视频二区| 成人在线观看不卡| 亚洲AV无码久久精品色欲| 精品久久人人爽人人玩人人妻| 亚洲国产成人精品无码区性色| 国产女人综合久久精品视| 国产精品无码一二三视频| 综合社区亚洲熟妇p| 2021国产精品自产拍在线| 丝袜国产一区| 免费精品一区二区h| 免费国产一级 片内射老| 91免费观看视频| 91成人精品视频| 91年精品国产福利线观看久久|