亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        考慮彈性影響的乘波體飛行動(dòng)力學(xué)特性

        2012-08-07 10:52:18向錦武曾開(kāi)春
        關(guān)鍵詞:配平馬赫數(shù)氣動(dòng)

        向錦武 曾開(kāi)春 聶 璐

        (北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院,北京100191)

        乘波體構(gòu)型近空間高超聲速飛行器(乘波體飛行器)可用于天地往返運(yùn)輸系統(tǒng)和高速導(dǎo)彈,在軍事和民用領(lǐng)域都具有廣闊的應(yīng)用前景.與常規(guī)飛行器相比,乘波體飛行器采用機(jī)體/發(fā)動(dòng)機(jī)一體化設(shè)計(jì),造成氣動(dòng)/推進(jìn)的強(qiáng)烈耦合[1].同時(shí),由于其細(xì)長(zhǎng)體外形的機(jī)體結(jié)構(gòu)中使用了大量輕質(zhì)材料,機(jī)體結(jié)構(gòu)柔性大、振動(dòng)頻率低,在氣動(dòng)載荷激勵(lì)下,可能發(fā)生顯著的彈性變形和彈性振動(dòng).機(jī)體結(jié)構(gòu)的大幅度彈性變形不僅改變了全機(jī)的氣動(dòng)特性,還引起推進(jìn)系統(tǒng)空氣流量以及入流參數(shù)的變化,形成乘波體飛行器特有的氣動(dòng)/彈性/推進(jìn)耦合現(xiàn)象[2-3].而機(jī)體的低頻彈性振動(dòng)則可能激發(fā)剛體動(dòng)力學(xué)短周期模態(tài),形成剛?cè)狁詈犀F(xiàn)象.因此,結(jié)構(gòu)彈性加劇了乘波體飛行動(dòng)力學(xué)的耦合性,并使其呈現(xiàn)出多尺度、強(qiáng)非線(xiàn)性的特點(diǎn),進(jìn)而對(duì)總體設(shè)計(jì)和控制系統(tǒng)設(shè)計(jì)提出了嚴(yán)峻挑戰(zhàn).由于結(jié)構(gòu)彈性參數(shù)會(huì)因設(shè)計(jì)更改而改變,并受到氣動(dòng)加熱等因素的影響,其值可能在較大范圍變化.因此,有必要研究結(jié)構(gòu)彈性參數(shù)對(duì)乘波體飛行器多耦合飛行動(dòng)力學(xué)特性的影響規(guī)律.

        乘波體飛行器動(dòng)力學(xué)特性分析的關(guān)鍵是建立多耦合動(dòng)力學(xué)模型.目前,主要的建模方法包括基于FEM(Finite Elenent Melthod)結(jié)構(gòu)模態(tài)計(jì)算的動(dòng)力學(xué)建模[4]、基于 CFD(Computationd Fluid Dynamics)數(shù)值擬合的動(dòng)力學(xué)建模[5]、基于近似理論方法的耦合動(dòng)力學(xué)建模等[6-7].前2種方法主要面向工程,其分析過(guò)程需要詳細(xì)的結(jié)構(gòu)布置與外形尺寸參數(shù),而近似理論建模方法主要面向科學(xué)研究,其計(jì)算效率相對(duì)較高,是目前乘波體飛行器飛行動(dòng)力學(xué)研究的主要建模方法.因此,本文采用近似理論建模方法,將氣動(dòng)外形、推進(jìn)系統(tǒng)、結(jié)構(gòu)模型等進(jìn)行了簡(jiǎn)化處理,再分別利用近似理論氣動(dòng)力估算方法、瑞利流理論、假設(shè)模態(tài)法等建立了考慮多種耦合因素的動(dòng)力學(xué)模型.研究了在不同機(jī)體結(jié)構(gòu)剛度下,結(jié)構(gòu)彈性對(duì)飛行動(dòng)力學(xué)靜配平特性及動(dòng)態(tài)特性的影響規(guī)律,為總體設(shè)計(jì)及飛行控制系統(tǒng)設(shè)計(jì)提供了參考.

        1 高超聲速飛行器飛行動(dòng)力學(xué)模型

        典型的乘波體飛行器采用機(jī)體/發(fā)動(dòng)機(jī)一體化設(shè)計(jì).機(jī)體前部呈尖楔形,高超聲速飛行時(shí),產(chǎn)生附體斜激波,為飛行器提供升力,同時(shí)為發(fā)動(dòng)機(jī)提供高壓空氣;機(jī)體后部是發(fā)動(dòng)機(jī)的外噴管,產(chǎn)生額外的推力,同時(shí)也產(chǎn)生附加的升力和俯仰力矩.因此,動(dòng)力學(xué)建模需要反映氣動(dòng)與推進(jìn)的這種耦合關(guān)系.本文分別采用激波膨脹波方法、瑞利流理論、假設(shè)模態(tài)法等對(duì)氣動(dòng)、推進(jìn)系統(tǒng)和結(jié)構(gòu)彈性動(dòng)力學(xué)進(jìn)行建模,給出了考慮氣動(dòng)/推進(jìn)/結(jié)構(gòu)耦合的乘波體飛行器動(dòng)力學(xué)模型.

        1.1 氣動(dòng)力模型

        采用基于激波膨脹波方法[7]的當(dāng)?shù)亓骰钊碚摚?]給出高超聲速非定常氣動(dòng)參數(shù),進(jìn)而建立了乘波體飛行器的氣動(dòng)模型.激波膨脹波方法適用于附體激波流動(dòng),能給出激波位置,適用的馬赫數(shù)范圍大,是計(jì)算高超聲速附體激波流動(dòng)氣動(dòng)力的基本方法.而牛頓沖擊流理論[6]形式簡(jiǎn)單,可用于快速估算操縱面的氣動(dòng)力,一般適用于馬赫數(shù)5以上.當(dāng)?shù)亓骰钊碚撃軌蚩紤]非定常效應(yīng),其應(yīng)用范圍為馬赫數(shù)1.2~10,但由于高超聲速流動(dòng)時(shí)的非定常效應(yīng)影響較小,其使用限制可適當(dāng)放寬.當(dāng)?shù)亓骰钊碚撊缡?1)、式(2)所示.

        其中,Vn為氣動(dòng)微元瞬時(shí)法向擾動(dòng)速度;q為俯仰角速度;?wd/?x 和分別為彈性變形轉(zhuǎn)角和速度;aSE,pSE,VLocal分別為激波膨脹波方法計(jì)算的當(dāng)?shù)匾羲?、壓?qiáng)和氣流速度模量;利用當(dāng)?shù)亓骰钊碚撔拚膲簭?qiáng)分布,可以積分得到各個(gè)氣動(dòng)面的氣動(dòng)力和力矩.此外,由于考慮了彈性變形的影響,以上建立的氣動(dòng)力模型能夠反映氣動(dòng)與結(jié)構(gòu)彈性的耦合.

        1.2 推進(jìn)系統(tǒng)模型

        根據(jù)文獻(xiàn)[6]的超燃沖壓發(fā)動(dòng)機(jī)模型,可以將發(fā)動(dòng)機(jī)分成擴(kuò)散段、燃燒室和內(nèi)噴管.在激波膨脹波方法給出上游流動(dòng)參數(shù)和激波位置的基礎(chǔ)上,可以計(jì)算出推進(jìn)系統(tǒng)的入流參數(shù).再利用瑞利增能流[9]理論可以給出內(nèi)流參數(shù),具體計(jì)算方法參看文獻(xiàn)[6].在分析內(nèi)流的基礎(chǔ)上,根據(jù)發(fā)動(dòng)機(jī)入流和出口參數(shù),可給出推力的公式如下:

        其中,下標(biāo)i表示入口流動(dòng)參數(shù);下標(biāo)e表示出口參數(shù);下標(biāo)∞表示自由流參數(shù);為發(fā)動(dòng)機(jī)氣體流量.可以發(fā)現(xiàn),發(fā)動(dòng)機(jī)的狀態(tài)強(qiáng)烈依賴(lài)于入流參數(shù),而入流參數(shù)又決定于上游機(jī)體前部下表面的流動(dòng)參數(shù)與激波位置,因此該建模方法能夠考慮氣動(dòng)與推進(jìn)的耦合.

        1.3 機(jī)體結(jié)構(gòu)彈性模型

        乘波體飛行器剛度分布中部最大,并向兩端遞減,因此,彈性變形主要集中在兩端.在缺乏數(shù)據(jù)的設(shè)計(jì)初期,這種彈性變形方式適合采用基于彈性集中假設(shè)的兩段模型[6].該模型在機(jī)身前部和后部連接的中部引入彈性鉸,并假設(shè)變形時(shí),前部和后部均剛性轉(zhuǎn)動(dòng),彈性鉸發(fā)生變形.此時(shí),機(jī)體彈性變形位移wd可表示為模態(tài)廣義坐標(biāo)與模態(tài)函數(shù)的乘積:

        其中,由于僅考慮一階模態(tài),因此i=1.根據(jù)拉格朗日方程可以得出質(zhì)量歸一化的彈性模態(tài)函數(shù)為

        其中,J1,J2分別為重心前后機(jī)身段的俯仰轉(zhuǎn)動(dòng)慣量;L1,L2分別為飛行器前后機(jī)身段的長(zhǎng)度.在模態(tài)坐標(biāo)下,引入彈性模態(tài)廣義力來(lái)計(jì)算結(jié)構(gòu)彈性變形,其表達(dá)式為

        其中,fz為機(jī)體法向分布載荷;Fzj為法向集中載荷.機(jī)體前后部分中線(xiàn)的等效傾角改變量Δ1,Δ2與模態(tài)廣義坐標(biāo)的關(guān)系為

        高超聲速飛行器巡航飛行時(shí),氣動(dòng)加熱對(duì)結(jié)構(gòu)力學(xué)特性影響較大,但相對(duì)于氣動(dòng)彈性響應(yīng),這種影響是一個(gè)緩慢的過(guò)程,因此本文將氣動(dòng)加熱的影響考慮成結(jié)構(gòu)剛度在一定范圍內(nèi)變化.

        2 縱向飛行力學(xué)特性分析方法

        從基本的拉格朗日方程出發(fā),推導(dǎo)了縱向剛?cè)狁詈系娘w行動(dòng)力學(xué)模型.機(jī)體坐標(biāo)系下,柔性飛行器縱向運(yùn)動(dòng)學(xué)和耦合動(dòng)力學(xué)方程如下:

        其中,狀態(tài)變量h為高度;u為x向飛行速度;w為z向速度;θ為俯仰角;q為俯仰角速度;η為彈性模態(tài)廣義坐標(biāo)為彈性模態(tài)廣義速度;下文分析中使用的迎角α可以近似由w與u的比值給出;Fx,F(xiàn)z分別為x向和z向合力;M為俯仰合力矩;N為彈性廣義力.

        2.1 靜態(tài)配平與飛行包線(xiàn)

        求解靜配平問(wèn)題時(shí),需要首先對(duì)方程式(8)~式(13)進(jìn)行處理,即令飛行速度分量的時(shí)間導(dǎo)數(shù)和、俯仰角速度 q及其對(duì)時(shí)間的導(dǎo)數(shù)、彈性變形速度均為0,將原來(lái)的6個(gè)方程化簡(jiǎn)為5個(gè).然后在指定飛行高度和速度,利用Newton-Raphson方法求解關(guān)于迎角α、俯仰角θ、升降舵偏轉(zhuǎn)角δe(向下偏轉(zhuǎn)為正)、發(fā)動(dòng)機(jī)油門(mén)δT以及彈性廣義坐標(biāo)η等5個(gè)變量的非線(xiàn)性方程組,可獲得靜平衡解.

        對(duì)于飛行包線(xiàn),考慮以下3種邊界條件:

        1)最大動(dòng)壓邊界:熱防護(hù)系統(tǒng)的耐受范圍決定了飛行的最大動(dòng)壓,其準(zhǔn)確值需要由試驗(yàn)才能給出,本文取為95.76 kPa(2000 lb/ft).

        2)發(fā)動(dòng)機(jī)熱阻塞邊界:當(dāng)發(fā)動(dòng)機(jī)燃燒室出口馬赫數(shù)低于1時(shí),發(fā)動(dòng)機(jī)會(huì)因熱阻塞而不能正常工作.該邊界是由可壓縮流體力學(xué)理論[9]推導(dǎo)得出的,本文給出一定的裕度,即要求發(fā)動(dòng)機(jī)燃燒室出口馬赫數(shù)不得低于1.3.

        3)最大油門(mén)邊界:由于進(jìn)入發(fā)動(dòng)機(jī)的空氣有限,因此當(dāng)其中的氧氣完全與燃油發(fā)生化學(xué)反應(yīng),即等效油門(mén)為1后,不能通過(guò)增大油門(mén)來(lái)增加推力.同時(shí)考慮到進(jìn)入發(fā)動(dòng)機(jī)的氧氣并不能都發(fā)生化學(xué)反應(yīng),因此等效油門(mén)最大值取為0.9.

        上述的3個(gè)邊界條件中,條件1)確定了可配平飛行的最大速度和最低高度,其與飛行器參數(shù)無(wú)關(guān),僅取決于大氣環(huán)境和飛行速度;條件2)和條件3)確定了最低飛行速度和最大高度,二者均與飛行器的參數(shù)相關(guān).本文在求解條件2)和條件3)時(shí),首先選取了用于分析求解的馬赫數(shù)值;然后分別針對(duì)每個(gè)馬赫數(shù)值,使用二分法求解最大可配平高度;最后在速度-高度圖中,將各個(gè)離散的邊界點(diǎn)連接起來(lái)形成飛行包線(xiàn).

        2.2 動(dòng)態(tài)特性

        在配平狀態(tài)下對(duì)方程式(8)~式(13)做線(xiàn)性化處理,可以得到如式(14)、式(15)所示的線(xiàn)性動(dòng)力學(xué)分析模型.

        3 機(jī)體剛度對(duì)飛行力學(xué)特性的影響

        根據(jù)文獻(xiàn)[7]給出的典型乘波體飛行器參數(shù),從飛行動(dòng)力學(xué)靜配平特性、飛行包線(xiàn)要求、動(dòng)態(tài)特性3個(gè)角度研究了機(jī)體結(jié)構(gòu)剛度對(duì)飛行動(dòng)力學(xué)特性的影響規(guī)律.

        3.1 機(jī)體結(jié)構(gòu)剛度對(duì)靜態(tài)特性的影響

        機(jī)體一階彎曲自然頻率ω1代表了機(jī)體彎曲剛度的大小,其值越大表明機(jī)體彎曲剛度越大.以下分析了配平參數(shù)隨ω1的變化關(guān)系,其中ω1為圓頻率,單位為rad/s.

        圖1 前體彈性變形等效傾角隨馬赫數(shù)的變化關(guān)系

        圖2 迎角配平值隨馬赫數(shù)的變化關(guān)系

        圖3 配平狀態(tài)下的全機(jī)升阻比隨馬赫數(shù)的變化關(guān)系

        圖4 油門(mén)配平值比隨馬赫數(shù)的變化關(guān)系

        圖5 升降舵配平值比隨馬赫數(shù)的變化關(guān)系

        圖1~圖5為H=25 908 m(85 000 ft),Ma=8,ω1分別為無(wú)窮大(剛體),25,18,14 rad/s時(shí),定常平飛狀態(tài)參數(shù)和控制變量隨馬赫數(shù)的變化.可以看出:

        1)在高馬赫數(shù)下,結(jié)構(gòu)彈性對(duì)飛行器的配平狀態(tài)有較大的影響.對(duì)于同一機(jī)體剛度條件下的柔性飛行器,隨著馬赫數(shù)的增大,結(jié)構(gòu)彈性變形引起的機(jī)體前部?jī)A角急劇增大,甚至與設(shè)計(jì)角度接近,因而對(duì)靜配平特性產(chǎn)生的影響不可忽略.

        2)機(jī)體剛度下降會(huì)導(dǎo)致配平迎角減小.從圖1、圖2可以看出,在高馬赫數(shù)段,隨著結(jié)構(gòu)剛度減小,結(jié)構(gòu)變形引起機(jī)身前部下表面傾角增大,進(jìn)而導(dǎo)致所需的配平迎角急劇減小.

        3)機(jī)體剛度下降會(huì)導(dǎo)致飛行性能降低.從圖3、圖4可以看出,隨著機(jī)體剛度下降,飛行器定常飛行時(shí)的升阻比將減小,所需的油門(mén)增大,并且在高馬赫數(shù)段,該影響變得更加顯著.

        4)機(jī)體剛度下降會(huì)導(dǎo)致升降舵配平規(guī)律改變.從圖5可以發(fā)現(xiàn),在所有馬赫數(shù)段,剛體飛行器的升降舵配平值隨馬赫數(shù)增大而減小,為正操縱;而對(duì)于ω1<25 rad/s的柔性飛行器,在低馬赫數(shù)區(qū)段,升降舵配平值隨馬赫數(shù)增大而增大,為正操縱;在高馬赫數(shù)區(qū)段,升降舵配平值隨馬赫數(shù)增大而減小,變?yōu)榉床倏v.因此控制系統(tǒng)設(shè)計(jì)時(shí),有必要檢驗(yàn)柔性高超聲速飛行器是否存在高馬赫數(shù)區(qū)反操縱現(xiàn)象.

        圖6 是 ω1分別為∞,18,14,10rad/s時(shí),乘波體飛行器的飛行包線(xiàn).可以看出機(jī)體剛度減小會(huì)顯著縮小高馬赫數(shù)段的可配平區(qū)域.隨著ω1的減小,發(fā)動(dòng)機(jī)熱阻塞邊界向馬赫數(shù)升高方向移動(dòng),而最大油門(mén)邊界向高度減小方向移動(dòng),隨著馬赫數(shù)升高,其減小量增大.此外,還可以發(fā)現(xiàn)在 ω1=10 rad/s時(shí),在馬赫數(shù)16附近,最大油門(mén)邊界存在不連續(xù)現(xiàn)象.分析其產(chǎn)生的機(jī)理首先要從高度對(duì)配平特性的影響入手,限于篇幅,本文直接給出結(jié)論:高度降低,剛體飛行器配平迎角減小,配平阻力減小,并且結(jié)構(gòu)彈性變形增大.在該結(jié)論的基礎(chǔ)上,可以發(fā)現(xiàn),由于結(jié)構(gòu)變形引起阻力增大,因此,當(dāng)彈性變形引起的阻力增加與配平迎角減小引起的阻力減小相抵消時(shí),通過(guò)降低高度的方式不能配平飛行器.

        圖6 不同機(jī)身彎曲剛度下的飛行包線(xiàn)

        3.2 機(jī)體剛度對(duì)動(dòng)態(tài)特性的影響

        在H=25908 m,Ma=8狀態(tài)下,分別計(jì)算ω1為∞,25,18,14 rad/s,升降舵為輸入,軌跡角為輸出的系統(tǒng)開(kāi)環(huán)的零極點(diǎn),如圖7所示.可以看出,高超聲速飛行器縱向飛行動(dòng)力學(xué)存在3個(gè)模態(tài),由5個(gè)特征根組成.其中短周期特征根由絕對(duì)值較大的一正實(shí)根和一負(fù)實(shí)根組成,表明其短周期缺少回復(fù)力,為靜發(fā)散.隨著機(jī)體彎曲剛度的減小,正實(shí)根進(jìn)一步增大,加劇了短周期的不穩(wěn)定.長(zhǎng)周期模態(tài)由2個(gè)小阻尼、低頻率共軛復(fù)根組成,隨著機(jī)體彎曲剛度減小,其實(shí)部變化較小,而虛部向?qū)嵼S靠近.高度模態(tài)對(duì)應(yīng)一個(gè)絕對(duì)值接近0的實(shí)根,其隨機(jī)體剛度變化較小.此外,還可以看出,機(jī)體彎曲剛度對(duì)零點(diǎn)位置的影響較小.綜上所述,機(jī)體剛度降低,短周期模態(tài)穩(wěn)定性急劇下降,進(jìn)而使按照剛體模型設(shè)計(jì)的控制系統(tǒng)失效.

        圖7 開(kāi)環(huán)零極點(diǎn)隨機(jī)身剛度減小的變化

        4 結(jié)論

        通過(guò)對(duì)考慮機(jī)體彈性的乘波體飛行器的縱向飛行動(dòng)力學(xué)建模、計(jì)算,分析了靜態(tài)配平特性、飛行包線(xiàn)及動(dòng)穩(wěn)定性隨機(jī)體彎曲剛度的變化規(guī)律,結(jié)果表明:

        1)隨著機(jī)體剛度減小,配平升阻比降低,耗油率增大,由發(fā)動(dòng)機(jī)熱阻塞和最大油門(mén)確定的飛行包線(xiàn)縮小;

        2)當(dāng)機(jī)體剛度較小時(shí),飛行器升降舵配平規(guī)律在不同馬赫數(shù)區(qū)段將發(fā)生變化,在低馬赫數(shù)段升降舵操縱為正操縱,在高馬赫數(shù)段變?yōu)榉床倏v;

        3)縱向存在不穩(wěn)定的短周期模態(tài),機(jī)體剛度減小,長(zhǎng)周期模態(tài)頻率將降低,短周期模態(tài)正實(shí)特征根增大,從而使系統(tǒng)穩(wěn)定性降低.

        References)

        [1]Fidan B,Mirmirani M,Ioannou P.Flight dynamics and control of air-breathing hypersonic vehicles:review and new Directions[R].AIAA-2003-708,2003

        [2]Cockrell C E,Engelund W C.Integrated aeropropulsive computational fluid dynamics methodology for the hyper-X flight experiment[J].Journal of Spacecraft and Rockets,2001,38(6):836 -843

        [3]黃喜元,王青,董朝陽(yáng).基于動(dòng)態(tài)逆的高超聲速飛行器魯棒自適應(yīng)控制[J].北京航空航天大學(xué)學(xué)報(bào),2011,37(5):560-563 Huang Xiyuan,Wang Qing,Dong Chaoyang.Robust adaptive control for hypersonic vehicle based on dynamic inversion[J].Journal of Beijing University of Aeronautics and Astronautics,2011,37(5):560 -563(in Chinese)

        [4]Heeg J,Zeiler T A,Pototzky AS,et al.Aerothermoelastic analysis of a NASP demonstrator model[R].NASA-TM-109007,1993

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

        [6]Chavez F R,Schmidt D K.Analytical Aeropropulsive/Aeroelastic hypersonic vehicle model with dynamic analysis[J].AIAA Journal of Guidance,Control,and Dynamics,1994,17(6):1308 -1319

        [7]Bolender M A,Doman D B.Nonlinear longitudinal dynamical model of an air-Breathing hypersonic vehicle[J].Journal of Spacecraft and Rockets,2007,44(2):374 -387

        [8]McNamara J J,Crowell A R,F(xiàn)riedmann P R,et al.Approximate modeling of unsteady aerodynamics for hypersonic aeroelasticity[J].Joural of Aircraft,2010,47(6):1932 -1945

        [9]Anderson J D.Modern compressible flow[M].2nd ed.New York:McGraw-Hill College,1995:77 -91

        [10]方振平,陳萬(wàn)春,張曙光.航空飛行器飛行動(dòng)力學(xué)[M].北京:北京航空航天大學(xué)出版社,2005:219-220 Fang Zhenping,Chen Wanchun,Zhang Shuguang.Aircraft flight dynamics[M].Beijing:Beihang University Press,2005:219 -220(in Chinese)

        猜你喜歡
        配平馬赫數(shù)氣動(dòng)
        高馬赫數(shù)激波作用下單模界面的Richtmyer-Meshkov不穩(wěn)定性數(shù)值模擬
        爆炸與沖擊(2024年7期)2024-11-01 00:00:00
        中寰氣動(dòng)執(zhí)行機(jī)構(gòu)
        一維非等熵可壓縮微極流體的低馬赫數(shù)極限
        基于NACA0030的波紋狀翼型氣動(dòng)特性探索
        載荷分布對(duì)可控?cái)U(kuò)散葉型性能的影響
        配平化學(xué)方程式小竅門(mén)——“單質(zhì)最后配平法”
        化學(xué)方程式的配平方法
        基于反饋線(xiàn)性化的RLV氣動(dòng)控制一體化設(shè)計(jì)
        化合價(jià)歸零法配平復(fù)雜氧化還原反應(yīng)方程式
        B737NG飛機(jī)安定面配平非典型故障分析
        人妻无码中文专区久久综合| 成人精品一区二区三区电影| 国产真实强被迫伦姧女在线观看| 视频一区视频二区制服丝袜| 日韩精品成人无码专区免费| 日本无遮挡吸乳呻吟视频| 日韩国产欧美视频| 日本精品一区二区在线看| 国内嫩模自拍偷拍视频| 国产人妻鲁鲁一区二区| 午夜精品一区二区三区的区别| 中国精学生妹品射精久久| 中文字幕大乳少妇| 国产三区二区一区久久| 在线观看成人无码中文av天堂| 亚洲色偷偷色噜噜狠狠99| 人妻丝袜中文字幕久久| 国产性感丝袜在线观看| 国产精品美女久久久久久| 亚洲精品有码在线观看| 92自拍视频爽啪在线观看| 精品亚洲麻豆1区2区3区| 色婷婷久久综合中文久久蜜桃av| 国产中文字幕乱码在线| 亚洲国产一区二区av| 国产精品无码久久综合| 欧美国产精品久久久乱码| 国产女人91精品嗷嗷嗷嗷| 日本免费播放一区二区| 欧美最猛黑人xxxx| 亚洲另类激情综合偷自拍图 | 丁香美女社区| 国产精品久久久久久久成人午夜| 青青草原亚洲在线视频| 自拍视频在线观看国产| 天堂中文官网在线| 黄视频国产| 国产一区二区三区不卡在线播放| 很黄很色很污18禁免费| 免费av片在线观看网站| 一亚洲一区二区中文字幕|