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

飛機起落架動力學建模及地面運動仿真

2021-07-12 04:38:58唐瑞琳
科學技術(shù)與工程 2021年16期
關(guān)鍵詞:飛機方向

唐瑞琳, 鞏 磊, 王 博

(航空工業(yè)第一飛機設(shè)計研究院飛控系統(tǒng)設(shè)計研究所, 西安 710089)

在飛機研制的過程中,需要在飛行模擬器上進行飛行品質(zhì)評估試驗以評價飛行品質(zhì),并在后期培訓飛行員以降低試飛風險。除空中飛行外,飛行模擬器還需要在飛機滑跑、離地、觸地階段能準確地模擬飛機的地面運動特性,提供給飛行員真實的起降感受。由于地面存在支反力和摩擦力,飛機在空中和地面的受力情況有著本質(zhì)區(qū)別,不僅導致飛行員在空中和地面的操縱方式不同,對于電傳飛機而言,一般也是依據(jù)起落架承載情況,切換飛控系統(tǒng)控制律的邏輯或結(jié)構(gòu)。因此對起落架動力學模型的要求不僅是起著地面支反作用,而是應盡可能地精準。特別對于螺旋槳飛機,由于螺旋槳同向旋轉(zhuǎn)引起不對稱滑流,地面滑跑時會受到滑流和氣動的側(cè)力,需要飛行員通過偏轉(zhuǎn)轉(zhuǎn)向輪(或方向舵)進行糾偏,因此需要正確描述起落架在橫側(cè)向的受力情況,才能反映飛機地面運動特性。這對起落架動力學模型精準度的要求更加嚴格。由于無人機自主滑跑仿真的需要,通過對無人機起落架的建模[1]可用于分析無人機特定的自主滑跑仿真[2-7]及轉(zhuǎn)彎操縱分析[8]。

由于起落架結(jié)構(gòu)、受力、工況復雜,通過建立虛擬樣機在ADAMS環(huán)境中仿真是一種行之有效的起落架系統(tǒng)設(shè)計手段[9-13],但飛行模擬器具有高實時性要求,因此虛擬樣機的方法不適用。在以往起落架實時仿真模型中,起落架結(jié)構(gòu)和受力的簡化和假設(shè)造成了一定的局限性。楊春生等[14]所建立的起落架動力學模型只限于特殊姿態(tài)的解算;Sun等[15]只在二維對稱平面建立了起落架動力學模型;蘇彬等[16]建模假設(shè)了起落架前輪轉(zhuǎn)向軸垂直安裝并忽略機輪傾側(cè)效應;徐冬苓等[17]和徐世玥[18]忽略了起落架的側(cè)向受力;張雯等[19]獨立計算機輪和緩沖器受力導致受力不平衡;袁東[20]未考慮俯仰姿態(tài)對機輪觸地判斷的影響且著陸時緩沖器壓縮量高頻振蕩。

為了建立能真實仿真飛機任意狀態(tài)的地面運動特性的起落架動力學模型以更好的評價飛行品質(zhì)和培訓飛行員,首先提出極坐標結(jié)合坐標轉(zhuǎn)換的方法計算每個起落架機輪的觸地點并得到承載情況,然后推導緩沖器壓縮行程和壓縮速率的代數(shù)計算方法,再分析得到機輪在空間的滾動方向定義,從而解析得到精準的側(cè)向摩擦,并利用相互作用力原理在不忽略緩沖器應力的基礎(chǔ)上解算出緩沖器受力及地面摩擦力,最后得到所有起落架相對于飛機的力和力矩,用于飛機運動參數(shù)解算,從而可模擬飛機任意狀態(tài)下的地面運動。

1 模型結(jié)構(gòu)

飛行模擬器飛行仿真模型的原理是將飛機的氣動、發(fā)動機和起落架在飛機體軸系下的力和力矩代入到體軸系下的運動方程進行統(tǒng)一解算,得到飛機運動參數(shù)。為了將起落架模型與仿真模型集成,起落架模型計算的結(jié)果應該是當前狀態(tài)下在飛機體軸系下的力和力矩,如圖1所示。

圖1 起落架動力學模型與飛行動力學模型的關(guān)系Fig.1 The relationship between landing gear dynamics model and flight dynamics model

飛機在地面運動狀態(tài)下,起落架緩沖器由于壓縮有行程和速度,形成了緩沖力,通過機輪作用到地面,地面再反作用給機輪支反力和摩擦力,從而形成起落架在體軸系的力和力矩。因此,起落架動力學模型首先需要根據(jù)飛機當前運動狀態(tài)計算起落架承載情況及緩沖器壓縮運動參數(shù),然后根據(jù)因壓縮形成的緩沖力解算機輪受到的支反力和摩擦力,最后將機輪受力投影到飛機體軸系上,輸出體軸系上的力和力矩供運動方程解算,其結(jié)構(gòu)示意圖如圖2所示,其中前輪偏轉(zhuǎn)角度、前輪輪軸與轉(zhuǎn)向軸的向量為前起落架所特有結(jié)構(gòu);最后將所有起落架在體軸系上的力和相對于重心的力矩進行疊加,得到起落架作用于飛機的力和力矩。

圖2 起落架動力學模型結(jié)構(gòu)Fig.2 The structure of landing gear dynamics model

為確保起落架動力學模型的準確性,需根據(jù)起落架原本結(jié)構(gòu)進行建模,起落架模型結(jié)構(gòu)如圖3所示。

圖3中,λ為前起落架前傾角,pm和pn分別為主起落架和前起落架在機體構(gòu)造坐標系下壓縮方向的向量,(xO,yO,zO)和(xo,yo,zo)分別為主起落架和前起落架未壓縮時輪軸在機體構(gòu)造系的坐標,R和r分別為主起落架和前起落架機輪的半徑,ln為前起落架輪軸到壓縮軸在機體構(gòu)造坐標系下的向量。

2 輪載的判定

由于機輪外輪廓是圓形,機輪距離地面的最低點(或機輪觸地點)與飛機的姿態(tài)相關(guān),其在飛機構(gòu)造系(以飛機鼻錐為原點,x軸沿機身軸線向后為正,y軸垂直于飛機對稱平面向右為正,z軸垂直于x-y平面向上為正)下的坐標不固定,因此首先需要確定機輪最低點在構(gòu)造系下的坐標,從而得到機輪軸的地面高度。

由于飛機在偏航時,起落架機輪跟隨飛機進行偏航,所以起落架機輪最低點(觸地點)只與飛機俯仰和滾轉(zhuǎn)姿態(tài)相關(guān)(轉(zhuǎn)向起落架還與轉(zhuǎn)向角和轉(zhuǎn)向軸前/后傾角有關(guān))。先以主起落架為例,起落架機輪相對于地面的姿態(tài)角定義同飛機歐拉角一樣,先俯仰后滾轉(zhuǎn)。

2.1 主起落架

建立機輪構(gòu)造坐標系,以機輪軸心為原點,x軸在機輪對稱面內(nèi)平行于機體構(gòu)造系x-y面向后為正,z軸在機輪對稱面內(nèi)垂直機體構(gòu)造系x-y面向上為正,y軸垂直于輪胎對稱面向右為正,如圖4所示。

機輪外輪廓可由極坐標方程表示為

(1)

式(1)中:下標w表示機輪構(gòu)造坐標系坐標;R為主起落架機輪半徑;ε為極坐標因子,為原點到輪廓任意一點與x軸的夾角。

建立機輪牽連地軸系,以機輪輪軸為原點,x軸為機輪構(gòu)造坐標系x軸在地面的投影向后為正,z軸為垂直于地面向上為正,y軸垂直于x-z平面向右為正,如圖5所示。

圖5 機輪牽連地軸系示意圖Fig.5 Schematic diagram of wheel implicated earth coordinate system

若機輪對稱面和機體對稱面平行,則機輪外輪廓坐標轉(zhuǎn)換到機輪牽連地軸系下為

(2)

式(2)中:下標mg表示主起落架機輪牽連地軸系坐標;Lgw為機輪外輪廓坐標系到機輪牽連地軸系的轉(zhuǎn)換矩陣;L′wg為機輪牽連地軸系到機輪外輪廓坐標系轉(zhuǎn)換矩陣的轉(zhuǎn)置。

Lwg=Lx(φ)Ly(θ)

(3)

式(3)中:Li(σ)為繞i軸旋轉(zhuǎn)角σ的轉(zhuǎn)換矩陣,i可以為x、y、z軸中任一個;σ可以為φ、θ中任一個,φ為滾轉(zhuǎn)角;θ為俯仰角。

(4)

機輪最低點(觸地點)的地軸z坐標有最小值,將z軸坐標對ε求導,可得

(5)

因此機輪最低點(觸地點)在機體構(gòu)造系的坐標(xm,ym,zm)為

(6)

式(6)中:(xO,yO,zO)為起落架支柱未壓縮時,主起落架機輪輪軸在機體構(gòu)造系的坐標。

由飛機的姿態(tài)可推出機輪最低點(觸地點)在地軸系下的坐標(xmG,ymG,zmG),可表示為

(7)

式(7)中:(xg,yg,zg)為飛機重心在構(gòu)造坐標系下的坐標;(xG,yG,zG)為飛機重心在地軸系下的坐標;ψ為偏航角。

如果求得zmG≥0(地軸系向下為正),則該主起落架機輪觸地有輪載。

2.2 前起落架

前起落架機輪最低點(觸地點)的計算與主起落架類似,不同的是前起落架機輪可偏轉(zhuǎn),而且由于前起落架有前傾角的原因,前起落架機輪偏轉(zhuǎn)的旋轉(zhuǎn)軸和體軸系的z軸不平行,故需先將機輪構(gòu)造坐標(定義同主起落架機輪)轉(zhuǎn)換到和體軸系平行,經(jīng)前機輪偏轉(zhuǎn)后再轉(zhuǎn)換回前起落架的機輪構(gòu)造坐標系。因此對于前起落架機輪,機輪外輪廓坐標轉(zhuǎn)換到機輪牽連地軸系的過渡矩陣為

Lwg=Ly(-λ)Lz(δ)Ly(λ)Lx(φ)Ly(θ)

(8)

式(8)中:δ為前輪偏度(機輪左偏為正);λ為前起落架前傾角(支柱前伸為正)。

ε=arctan{(cosθcosφcosλ-sinθsinλ)cosλ+[(sinθcosλ+cosθcosφsinλ)cosδ-

cosθsinφsinδ]sinλ}{(cosθcosφcosλ-

sinθsinλ)sinλ-[(sinθcosλ+

cosθcosφsinλ)cosδ-cosθsinφsinδ]cosλ}-1

(9)

由于前起落架機輪輪軸一般不在轉(zhuǎn)向軸上,所以前輪偏轉(zhuǎn)時輪軸的機體構(gòu)造系坐標會發(fā)生變化,令輪軸到轉(zhuǎn)向軸連接處在機體構(gòu)造坐標系下的向量為ln(前起落架未偏轉(zhuǎn)時),則前起落架機輪最低點(觸地點)在機體構(gòu)造系的坐標(xn,yn,zn)為

(10)

式(10)中:(xo,yo,zo)為前起落架未轉(zhuǎn)向的機輪輪軸在起落架支柱未壓縮時的機體構(gòu)造系坐標。

再由飛機姿態(tài)推出機輪最低點(觸地點)在地軸系下的坐標(xnG,ynG,znG)可表示為

(11)

如果求得znG≥0(地軸系向下為正),則前起落架機輪觸地,有輪載。

3 起落架緩沖力的計算

起落架緩沖器所提供的緩沖力是起落架提供給飛機支反力的主要組成部分,主要由隨壓縮行程變化的彈簧力、隨壓縮速率變化的阻尼力以及緩沖器內(nèi)部有相對運動而產(chǎn)生的滑動摩擦力組成。

3.1 壓縮行程

當計算得到主起落架機輪觸地有輪載時,則主起落架緩沖器被壓縮。令主起落架的壓縮方向在機體構(gòu)造系下為pm=[pmx,pmy,pmz],將其轉(zhuǎn)換到地軸系下,壓縮方向pmG,可表示為

(12)

由此可得主起落架壓縮行程Lm為

(13)

上述計算過程對于前起落架同樣適用。

3.2 壓縮速率

壓縮速率可通過緩沖器無壓縮時(假設(shè)未觸地)機輪輪軸地速和實際壓縮時(實際已觸地)的機輪輪軸地速在壓縮方向的差求得。輪軸的地速是重心速度與繞重心轉(zhuǎn)動的線速度之和。

3.2.1 無壓縮時機輪輪軸沿壓縮方向速度

令飛機重心的地速為Vg=[Vgx,Vgy,Vgz],則無壓縮時主起落架機輪輪軸沿壓縮方向的地速Vmg_uncomp可表示為

(14)

3.2.2 壓縮時機輪輪軸沿壓縮方向速度

壓縮時機輪輪軸沿壓縮方向的速度計算方法與未壓縮時相似,不同點在于:①輪軸與重心在機體構(gòu)造系下的相對位置有變化;②壓縮時輪軸地速在機輪牽連地軸系z軸的分量為0。由此可推出主起落架壓縮時機輪輪軸沿壓縮方向的地速Vmg_comp可表示為

(15)

因此主起落架緩沖器的壓縮速率為

Vm_comp=Vmg_uncomp-Vmg_comp

(16)

同理可計算前起落架的壓縮速率,不同的是式(14)、式(15)中未壓縮機輪輪軸在機體坐標系的坐標(xO,yO,zO)修正為前起落架偏轉(zhuǎn)后的未壓縮機輪輪軸在機體坐標系的坐標,即

(17)

3.3 緩沖力

根據(jù)起落架緩沖器的結(jié)構(gòu)和參數(shù)可以得到緩沖器彈簧力、阻尼力與壓縮行程、壓縮速率的函數(shù)關(guān)系;滑動摩擦和壓縮速率有關(guān),即有壓縮速率就有反向的固定滑動摩擦力(滑動摩擦力具有消除仿真振蕩的作用,只采用彈簧-阻尼特性的緩沖器會存在振蕩現(xiàn)象[21])。由此可得主起落架的緩沖力Fm_comp為

Fm_comp=Fm_spring(Lm)-sgn(Vm_comp)×

Fm_damp(Lm,Vm_comp)-sgn(Vm_comp)Fm_0

(18)

式(18)中:Fm_spring(Lm)為主起落架緩沖器彈簧力與壓縮行程的函數(shù)關(guān)系;Fm_damp(Lm,Vm_comp)為阻尼力與壓縮行程和壓縮速率的函數(shù)關(guān)系;Fm_0為滑動摩擦由起落架緩沖器的設(shè)計決定,其具體形式由起落架設(shè)計方提供。同理可得前起落架的緩沖力。

4 機輪摩擦力的構(gòu)成

起落架機輪和地面的受力關(guān)系除了地面提供支反力以外還有摩擦力,摩擦力根據(jù)機輪滾動的方向分為縱向摩擦力和側(cè)向摩擦力,其中沿機輪滾動方向的摩擦分量為縱向摩擦,垂直于滾動方向的摩擦分量為側(cè)向摩擦。

4.1 縱向摩擦

主起落架的縱向摩擦包含滾動摩擦和剎車壓力摩擦,前起落架只有滾動摩擦。滾動摩擦系數(shù)根據(jù)跑道面的不同而為不同的常數(shù),剎車壓力摩擦系數(shù)是剎車壓力的函數(shù),由于剎車壓力與剎車踏板偏度有近似一一對應的關(guān)系,剎車壓力摩擦系數(shù)也可以認為是剎車踏板偏度的函數(shù)。由此,起落架的縱向摩擦系數(shù)為

(19)

式(19)中:μm_long為主起落架縱向摩擦系數(shù);μn_long為前起落架縱向摩擦系數(shù);μ0為滾動摩擦系數(shù);μbreke(pedal)為剎車摩擦系數(shù),其中pedal為剎車程度。

4.2 側(cè)向摩擦

當機輪運動方向不在機輪對稱平面內(nèi)時,機輪存在側(cè)向摩擦。機輪的側(cè)向摩擦系數(shù)是機輪運動方向與機輪滾動方向夾角(即側(cè)偏角βm)的函數(shù),側(cè)偏角較小時,側(cè)向摩擦系數(shù)與側(cè)偏角近似成線性關(guān)系。側(cè)偏角可視為機輪運動方向和機輪地面坐標系x軸的夾角與機輪地面坐標系x軸和機輪滾動方向夾角的和。此外由于機輪存在最大摩擦系數(shù),側(cè)向摩擦系數(shù)應小于最大摩擦系數(shù)與縱向摩擦系數(shù)的幾何差,所以主起落架機輪側(cè)向摩擦系數(shù)的形式為

(20)

(21)

式中:μm_lat為主起落架側(cè)向摩擦系數(shù);μm_max為主起落架最大摩擦系數(shù);μβ(βm)為側(cè)向摩擦隨側(cè)偏角的函數(shù);(Vmgx,Vmgy,Vmgz)為主起落架輪軸在地軸系的地速;δmg為主起落架機輪滾動方向與機輪牽連地軸系的夾角。

μβ(βm)有多種輪胎模型可使用,可采用魔術(shù)公式輪胎模型或Dugoff輪胎模型[22]。為求解主起落架機輪滾動方向相對于機輪地軸牽連系的偏度δmg,需明確機輪滾動方向。由于機輪在只受摩擦力的情況下,在機輪輪廓平面內(nèi)滾動,因此定義一個垂直于輪軸到觸地點向量并在機輪對稱平面內(nèi)的單位向量,其機輪構(gòu)造坐標系的坐標為[-sinε,0,cosε],該單位向量在機輪牽連地軸系的投影與機輪牽連地軸系x軸的夾角即為δmg,可表示為

(22)

式(22)中:(δmgx,δmgy,δmgz)為主起落架單位向量在機輪牽連地軸系的方向。

前起落架機輪側(cè)向摩擦系數(shù)與主起落架相似,但由于前起落架可偏轉(zhuǎn)的原因,式(23)中未壓縮機輪輪軸在機體坐標系的坐標修正為前起落架偏轉(zhuǎn)后的未壓縮機輪輪軸在機體坐標系的坐標為

(23)

式(23)中:βn為前起落架機輪的側(cè)偏角;(Vngx,Vngy,Vngz)為前起落架輪軸在地軸系的地速;δng為前起落架機輪滾動方向與機輪牽連地軸系的夾角;Ln為前起落架的壓縮行程。

前起落架機輪滾動方向相對于機輪地面坐標系偏度的計算中,也要考慮前起落架前傾角和轉(zhuǎn)彎偏度的因素,可表示為

(24)

式(24)中:(δngx,δngy,δngz)為前起落架單位向量在機輪牽連地軸系的方向。

5 支反力的求解

對于機輪觸地點,地面對于機輪的支反摩擦力和機輪施加于地面的壓力是一對大小相等、方向相反的力。根據(jù)上述章節(jié)已求出緩沖力和地面各摩擦力的系數(shù),其中緩沖力只是這對力在緩沖器壓縮方向的分力,垂直于壓縮方向還存在應力;地面各摩擦系數(shù)和地面壓力的關(guān)系可確定合力的方向,因此已知分力大小和合力方向,可求解出合力的大小和方向,以主起落架為例,合力大小Fm為

(25)

建立機輪地面坐標系,以機輪觸地點為原點,x軸為輪軸滾動方向向前為正,y軸為輪軸在地面的投影向右為正,z軸為垂直于地面向下為正,如圖6所示。可得到合力在主起落架機輪地面坐標系的各分力Fmg為

圖6 機輪地面坐標系示意圖Fig.6 Wheel earth coordinate sytem

(26)

前起落架機輪由于有轉(zhuǎn)向,式(25)中Lwg的形式[式(8)]與主起落架有區(qū)別。

6 起落架相對機體的力和力矩

將機輪地面坐標系的分力轉(zhuǎn)換到體軸系上,可得主起落架對機體的力FGear_m為

FGear_m=Lx(φ)Ly(θ)Fmg

(27)

前起落架有同樣的形式。

主起落架對機體的力矩MGear_m由主起落架對機體的力叉乘機輪相對重心在體軸系的坐標可得

(28)

起落架的力矩計算需要將式中未壓縮機輪輪軸在機體坐標系的坐標(xo,yo,zo)修正為前起落架偏轉(zhuǎn)后的未壓縮機輪輪軸在機體坐標系的坐標,如式(17)所示。

將所有起落架的力和力矩相加就得到起落架系統(tǒng)對飛機在體軸系上的力FGear和力矩MGear為

(29)

7 仿真驗證

起落架動力學模型主要實現(xiàn)飛機在地面的落震、轉(zhuǎn)彎、剎車功能,為驗證該起落架動力學模型的正確性,對以上3個功能進行仿真驗證。同時為了適用于有運動平臺的飛行模擬器,確保著陸時的過載與實際相符,對著陸工況進行仿真驗證。

7.1 落震

設(shè)置飛機初始高度稍高于地面,俯仰角為3°,呈稍抬頭的姿態(tài),確保主起落架先觸地;速度為0.5 m/s,避免迎角側(cè)滑角解算出現(xiàn)奇點。落震過程重心高度和俯仰角變化曲線如圖7所示。

由圖7可見,飛機在落震過程中,重心高度經(jīng)過三次振蕩后穩(wěn)定(由于起落架有高度,靜止后飛機重心高度為正值),并且由于抬頭的初始姿態(tài),飛機初始有俯仰的振蕩,最終快速收斂穩(wěn)定。這與實際試驗現(xiàn)象相符[23]。

7.2 轉(zhuǎn)彎

設(shè)置飛機在第10秒左蹬腳蹬,第14秒時收回,偏航角、滾轉(zhuǎn)角、前起落架機輪側(cè)偏角及3個起落架機輪的軌跡如圖8所示。

機輪軌跡是由空中向地面俯瞰圖8 轉(zhuǎn)彎過程中仿真結(jié)果Fig.8 Simulation result of turning process

由圖8可以看到,左蹬腳蹬后飛機航向向左偏轉(zhuǎn),腳蹬回中后航向保持在松腳蹬的時刻;轉(zhuǎn)彎過程中,飛機由于向右的離心力會產(chǎn)生輕微的右滾轉(zhuǎn);前機輪偏轉(zhuǎn)后,機輪速度方向很快靠近滾動方向。

7.3 剎車

設(shè)置初始速度為40 m/s,第5 s開始持續(xù)滿剎車,速度、俯仰角變化如圖9所示,開始剎車后飛機低頭、速度急速減小;因為機輪位于重心以下,剎車帶來低頭力矩。

圖9 剎車過程的仿真結(jié)果Fig.9 Simulation result of breaking process

7.4 著陸

設(shè)置初始速度為著陸速度,飛機略高于地面并以抬頭姿態(tài)進行著陸。著陸時重心高度、俯仰角及重心法向過載如圖10所示,著陸觸地后,重心高度和俯仰角經(jīng)過一次振蕩就快速收斂。由于仿真初始設(shè)置為“硬著陸”,所以最大法向過載約3.5g(g為重力加速度),該過載值允許飛行模擬器帶上運動平臺給飛行員提供著陸觸地時的過載感受。

圖10 著陸過程的仿真結(jié)果Fig.10 Simulation result of landing process

8 結(jié)論

(1)在沒有忽略起落架結(jié)構(gòu)尺寸和受力方式的基礎(chǔ)上建立了受力特性完整的起落架動力學模型。

(2)由仿真結(jié)果可知,飛機落震經(jīng)過3次緩沖振蕩后收斂穩(wěn)定;飛機的航向跟隨前機輪的偏轉(zhuǎn)而改變,前機輪回中后保持回中時刻的航向,轉(zhuǎn)彎過程中向轉(zhuǎn)彎相反的方向傾斜;剎車時飛機速度快速較小并產(chǎn)生低頭的姿態(tài)。

(3)著陸仿真的過載仿真允許將該模型加載到帶運動平臺的飛行模擬器中,使飛行模擬器可模擬任意狀態(tài)下的地面運動并提供給飛行員真實的起降滑跑體驗。

猜你喜歡
飛機方向
鷹醬想要“小飛機”
飛機失蹤
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
國航引進第二架ARJ21飛機
“拼座飛機”迎風飛揚
當代陜西(2019年11期)2019-06-24 03:40:28
乘坐飛機
主站蜘蛛池模板: 欧美精品成人| 精品久久久久久成人AV| 欧美成一级| 原味小视频在线www国产| 日本高清有码人妻| 久久婷婷五月综合色一区二区| 97色婷婷成人综合在线观看| 欧美色丁香| 亚洲水蜜桃久久综合网站| 人妻无码中文字幕第一区| 嫩草在线视频| 伊人色天堂| 亚洲国产一区在线观看| 日韩不卡免费视频| 成人精品亚洲| 成人毛片免费观看| 极品国产在线| 一本色道久久88综合日韩精品| 天天操精品| 99精品这里只有精品高清视频| 22sihu国产精品视频影视资讯| 欧美曰批视频免费播放免费| 精品午夜国产福利观看| 成年免费在线观看| 久久免费视频6| 亚洲免费毛片| 亚洲精品麻豆| 亚洲欧美激情小说另类| 久久亚洲日本不卡一区二区| 成人日韩欧美| 福利在线不卡| 国产91精品调教在线播放| 日韩欧美国产精品| 国产午夜精品一区二区三| 亚洲欧洲国产成人综合不卡| 亚洲综合狠狠| 熟女视频91| 欧美第二区| 国内精品伊人久久久久7777人| 国产在线观看99| 国产91熟女高潮一区二区| 色哟哟国产成人精品| jizz在线免费播放| 亚洲第一成人在线| 亚洲精品欧美重口| 538国产视频| 免费人成网站在线观看欧美| 成人欧美日韩| 日韩精品高清自在线| 亚洲色偷偷偷鲁综合| 蜜臀av性久久久久蜜臀aⅴ麻豆| 伊人91视频| 国产精品手机在线播放| 91在线高清视频| 亚洲乱码视频| 99re经典视频在线| 久久熟女AV| 91在线精品免费免费播放| 欧美成人手机在线观看网址| 国产色婷婷| 国产波多野结衣中文在线播放| 精品久久久久久久久久久| 三上悠亚一区二区| 毛片免费观看视频| 全免费a级毛片免费看不卡| 亚洲最猛黑人xxxx黑人猛交| 好吊色国产欧美日韩免费观看| 国产在线精品99一区不卡| lhav亚洲精品| 51国产偷自视频区视频手机观看| 啪啪永久免费av| 亚洲高清中文字幕在线看不卡| 99青青青精品视频在线| 久久国产免费观看| 视频二区中文无码| 日韩专区第一页| 国产精品无码AV中文| 色成人亚洲| 亚洲欧洲一区二区三区| 九九九精品成人免费视频7| 亚洲第一区精品日韩在线播放| 视频一本大道香蕉久在线播放 |