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

基于特征根有界攝動分析的高超聲速飛行器雙通道控制

2022-11-09 09:58:12楊雨宸張增輝閆佳寧張晶楊凌宇
北京航空航天大學學報 2022年10期
關鍵詞:模態(tài)

楊雨宸 張增輝 閆佳寧 張晶 楊凌宇

(1. 北京航空航天大學 自動化科學與電氣工程學院, 北京 100083; 2. 北京機電工程總體設計部, 北京 100854)

高超聲速飛行器飛行狀態(tài)變化劇烈,包線范圍跨度極大,是近年來各國爭相發(fā)展的前沿技術[1-2]。 在多種氣動控制方案中,以HTV-2 為代表的高超聲速飛行器采用了一種雙通道控制方案[3],其僅用一對升降副翼完成了飛行器的姿態(tài)控制,優(yōu)點是可以有效降低飛行器結構復雜度、減輕質(zhì)量;但對飛行器本身的氣動特性和控制方法提出了更高的要求,尤其對于具有強不確定性的高超聲速飛行器,該問題更為復雜。

與傳統(tǒng)飛控系統(tǒng)設計不同,高超聲速飛行器的高馬赫飛行速度使得其對各類參數(shù)變化非常敏感,也由此帶來了強不確定性、強耦合等問題,給其控制系統(tǒng)的設計帶來了巨大挑戰(zhàn)。 近年來,學者們進行了大量關于高超聲速飛行器姿態(tài)控制的研究,包括自適應、魯棒H∞、滑模、容錯和智能等多種控制方法也都被運用在了其控制系統(tǒng)設計上[4-7]。 而對于類HTV-2 氣動布局飛行器的雙通道控制問題,相關研究相對較少且主要集中在理論研究層面。 文獻[8]研究了非最小相位欠驅(qū)動高超聲速飛行器的二階動態(tài)滑??刂茊栴},通過系統(tǒng)分解的方法設計了滑??刂破?。文獻[9]從工程應用的角度出發(fā),針對氣動控制舵面僅有體襟翼的情況提出了一種基于欠驅(qū)動構型飛行器耦合特性的橫側向控制策略。 文獻[10]針對只有2 個舵的欠驅(qū)動再入飛行器,分別針對快回路和慢回路設計了超扭曲滑模控制器和分層滑??刂破?組成了具有抗飽和功能的姿態(tài)跟蹤控制器。 文獻[11]提出了基于輸出重定義的非最小相位動態(tài)逆控制方法,利用根軌跡法確定控制輸出組合系數(shù),實現(xiàn)系統(tǒng)零動態(tài)的分配。

可以看出,先進控制方法在穩(wěn)定性理論和仿真驗證方面取得了重要的進展,但仍未在工程設計界廣泛應用,如何將先進控制理論與工程需求相結合仍是一個重要的問題。

在飛行器控制工程設計中,基于極點和模態(tài)的概念仍然是工程上分析和設計的重要手段,將欠驅(qū)動問題與模態(tài)設計結合起來,并提供適用的魯棒性分析方法,對改善現(xiàn)有飛行控制設計流程和推進現(xiàn)代控制理論的工程化有較強的研究意義。

針對上述問題,本文以雙通道高超聲速飛行器為對象,以系統(tǒng)的閉環(huán)極點為結合點,提出基于特征根有界攝動分析的魯棒性評價手段,并以六自由度非線性對象為例進行了對比分析和驗證。首先,建立了NASA Winged-Cone 飛行器六自由度非線性模型;然后,給出了荷蘭滾模態(tài)與氣動系數(shù)之間的近似關系,并提出了2 種可行的雙通道控制方案;同時,提出了特征根靈敏度矩陣與特征根有界攝動系數(shù)的概念及魯棒性分析方法;最后,基于非線性模型進行了對比分析與驗證。

1 高超聲速飛行器非線性建模

本文所建立的高超聲速飛行器再入模型采用NASA 蘭利研究中心提供的公開高超聲速飛行器Winged-Cone 的氣動數(shù)據(jù)及總體參數(shù)[12],其氣動布局如圖1 所示。

圖1 Winged-Cone 氣動布局[12]Fig.1 Winged-Cone pneumatic layout[12]

Winged-Cone 采用三角機翼布局,尾部有左右獨立工作的升降副翼和帶有方向舵的垂直尾翼,前端還有可收縮的水平鴨翼。

在飛機的舵面使用上,由于本文主要針對無動力滑翔階段,水平鴨翼收入機體不參與控制;另外,為開展雙通道控制研究,將方向舵始終固定于平衡位置,飛機僅采用左右升降副翼完成姿態(tài)控制。

建模中所用到的其他假設如下:

1) 飛行器為理想剛體。

2) 無動力滑翔過程中質(zhì)量不變,且推進系統(tǒng)不工作,即推力為0。

3) 飛行器為面對稱,故近似認為慣性積Ixy和Iyz為0。

4) 地球為平面大地。

基于上述假設給出被控對象質(zhì)心運動學模型為[13]

式中:[Vχγxyz]T分別為速度、航跡方位角、航跡傾角,以及質(zhì)心在地球坐標系下的三軸位移分量;L、D、Y分別為升力、阻力、側力;m為飛行器質(zhì)量;g為當?shù)刂亓铀俣取?/p>

姿態(tài)動力學模型為

式中:[αβμpqr]T分別為迎角、側滑角、傾側角、滾轉角速率、俯仰角速率和偏航角速率;Ix、Iy、Iz為轉動慣量;Ixz為慣性積;LA、MA、NA分別為飛行器總的滾轉力矩、俯仰力矩、偏航力矩。

非線性的氣動力和氣動力矩計算分別如下:

式中:CX、CY和CZ為氣流系三軸力系數(shù);Cl、Cm和Cn為三軸力矩系數(shù),具體數(shù)值詳見文獻[12];Q為動壓;S為機翼參考面積;b為翼展;c為氣動弦長;Xcg為重心與取矩中心距離。

各參數(shù)具體數(shù)值如表1 所示。

表1 Winged-Cone 主要參數(shù)Table 1 Key parameters of Winged-Cone

方程(1) ~(4)共同組成了所用的高超聲速飛行器再入段的六自由度非線性模型。

2 荷蘭滾模態(tài)分析及反饋控制方案

雙通道控制的一個重要問題在于飛機的橫側向處于欠驅(qū)動狀態(tài),無法獨立完成橫向和航向的解耦控制,因此,如何利用單個控制量抑制荷蘭滾模態(tài)的影響,并完成對橫滾運動的跟蹤控制,成為雙通道控制的一個重要問題。

本節(jié)先從Jacobi 線性化得到的橫側向運動模型出發(fā),推導荷蘭滾模態(tài)的極點與氣動參數(shù)關系;再結合工程要求,確定合理的反饋量與反饋形式,進而給出改善橫側向閉環(huán)特性的可行方案。

2.1 橫側向模態(tài)特性分析

選取橫側向狀態(tài)對模型進行Jacobi 線性化,可以得到如下線性模型:

其中:Y′*、N′*和L′*分別為側力、偏航力矩、滾轉力矩對下角標變量的偏導數(shù);γ0和α0分別為當前配平工作點的航跡傾角和迎角;δa為等效副翼。

矩陣AL中各元素為

式中:v0、α0分別為線性化工作點處的空速、迎角;Cnβ為偏航力矩系數(shù)關于側滑角偏導數(shù);CYβ為側力關于側滑角偏導數(shù);Clβ為滾轉力矩系數(shù)對側滑角偏導數(shù);Cnr、Clp分別為偏航、滾轉力矩系數(shù)的偏導數(shù);Cnp、Clr為交叉力矩系數(shù)的偏導數(shù)。

對于面對稱飛行器,矩陣AL特征方程一般為

式中:ξ和ω分別為荷蘭滾模態(tài)阻尼比和固有頻率; -λr為滾轉阻尼模態(tài)極點; -λs為對應螺旋模態(tài)極點;In為n×n單位矩陣。

考慮到高超聲速飛行器在整個飛行包線內(nèi)氣動特性變化顯著,因此,選取滑翔段設計包線為α∈[ -2°,14°],馬赫數(shù)Ma∈[4,12],對模型的橫側向特性進行分析,結果如圖2 ~圖4 所示。

圖4 氣動系數(shù)隨馬赫數(shù)變化曲線Fig.4 Aerodynamic coefficients with Mach number

圖2 為迎角12°時螺旋模態(tài)和滾轉阻尼模態(tài)極點隨Ma變化的曲線。 可以看出,Winged-Cone飛機與傳統(tǒng)飛機類似,其螺旋模態(tài)對應非常小的實根,滾轉阻尼模態(tài)對應較大的負實根,極點值在[ -0.03, -0.035]之間,大約為螺旋模態(tài)極點模值的5 ~10 倍。

圖2 螺旋與滾轉阻尼模態(tài)極點Fig.2 Spiral mode and roll-damping mode poles

圖3 為荷蘭滾阻尼比隨迎角與馬赫數(shù)變化曲線。 可以看出,在整個包線內(nèi)荷蘭滾模態(tài)阻尼很小,約在[0.005,0.025]之間,必須采用控制方法加以改善。

圖3 荷蘭滾阻尼比隨迎角和馬赫數(shù)變化曲線Fig.3 Dutch roll damping ratio with angle of attack and Mach number

另外,為了分析矩陣AL中各元素對閉環(huán)特性的影響,給出了矩陣AL主要氣動參數(shù)隨馬赫數(shù)變化的曲線,如圖4 所示。 可以看出,在所有氣動系數(shù)中,N′β與L′β比其他元素大1 ~2 個數(shù)量級。

根據(jù)上述分析,可以發(fā)現(xiàn)此類飛行器具有以下特點:

1) 滾轉阻尼模態(tài)極點數(shù)值遠遠大于螺旋模態(tài)極點,即λs?λr。

2) 荷蘭滾阻尼比ξ?1,需要改善。

3)N′β與L′β比其他氣動參數(shù)大1 ~2 個數(shù)量級。4) 副翼對航向控制效率近似為0,即N′δa?0。

2.2 荷蘭滾模態(tài)近似特性計算方法

傳統(tǒng)飛行控制中對荷蘭滾模態(tài)的近似分析方法基于克萊姆法則,在簡化過程中忽略了交叉倒數(shù)項,如L′β、N′p等[14],對包含方向舵的飛行器可得到簡單可行的控制方案;但對雙通道控制飛機,由于缺乏方向舵改善荷蘭滾阻尼,必須找到通過副翼改善荷蘭滾阻尼的手段,本節(jié)采用多項式分解的方式,對多項式的每一系數(shù)分別討論。

橫側向特征方程(7)可展開表示為

相對計算誤差定義如式(21)所示。 為了驗證式(21)的準確性,引入文獻[14]中的傳統(tǒng)方法進行對比,對設計包線內(nèi)的荷蘭滾模態(tài)的ξ和ω進行了計算,比較結果如圖5 和圖6 所示。

圖5 飛行包線內(nèi)ξ 相對估計誤差Fig.5 Relative estimation error of ξ in flight envelope

圖6 飛行包線內(nèi)ω 相對估計誤差Fig.6 Relative estimation error of ω in flight envelope

式中:ξ^ 和^ω分別為ξ和ω的真值。

圖5 和圖6 中,紅線為各馬赫數(shù)下本文方法得到的相對誤差,藍線為傳統(tǒng)方法計算得到的相對誤差[14]。 從圖中誤差結果來看,在整個設計包線內(nèi),所推導的荷蘭滾阻尼比相對誤差都在21%以下,固有頻率相對誤差小于1%,估計精度明顯優(yōu)于文獻[14]中的荷蘭滾預測方法,這表明推導的荷蘭滾模態(tài)關系足以滿足研究的需要。

2.3 橫側向反饋控制方案分析

雙通道控制的難點在于缺少產(chǎn)生航向力矩的舵面,因此其荷蘭滾模態(tài)的改善不能采用常規(guī)方式。 本節(jié)根據(jù)荷蘭滾模態(tài)與氣動系數(shù)之間的關系來研究可行的反饋控制方案。

由式(20),并注意Y′β、N′r和L′β均為負,可得出以下結論:

1) 減小Y′β、N′r和N′β均能改善荷蘭滾阻尼,這也是常規(guī)飛行器改善荷蘭滾特性的主要手段,但對于雙通道控制飛行器,其副翼對偏航和側力的操縱效率極低,需要較大增益才能對這3 項產(chǎn)生明顯影響,考慮到舵偏限制,該方式在工程上不具備可行性。

2) 觀察式(20)實部中L′p相關項,當α0很小時,2 項近似相消,因此增大L′p對荷蘭滾模態(tài)的改善效率很低,但增大L′p對改善滾轉阻尼模態(tài)具有較強的效益[14]。

3) 增大L′β絕對值能改善荷蘭滾阻尼,且由于副翼對滾轉操縱效率較高,不存在舵偏過限的問題,因此工程上是可以實現(xiàn)的。

綜合上述結論,較為合理的反饋狀態(tài)量選取應為側滑角與滾轉角速率,給出“極點配置方案”,該方案的主要思路為:通過合適的增益反饋改變L′β與L′p的值,進而改善閉環(huán)系統(tǒng)和荷蘭滾特性。 按照前述反饋策略,配置反饋增益:

除了上述“極點配置方案”,工程上經(jīng)常采用一種通道隔離的解耦設計方法,即將被控狀態(tài)與動態(tài)性能差的狀態(tài)解耦,以降低影響。 對于滑翔段,其主要控制量為航跡傾角,若將其與荷蘭滾模態(tài)解耦,也可以達到設計要求,因此本文提出了“解耦控制方案”,其基本思想為:偏航通道利用對象本身靜穩(wěn)定性實現(xiàn)鎮(zhèn)定,不再追求荷蘭滾模態(tài)快速收斂;反饋目的變?yōu)橄酵ǖ缹L轉通道的影響,也就是通過反饋來降低L′β;對滾轉通道實現(xiàn)合理增穩(wěn)。

解耦控制方案的控制結構與極點配置方案相同,但K的設計過程需要滿足:

3 特征根有界攝動分析方法

對于高超聲速飛行器,其強不確定性會直接影響系統(tǒng)閉環(huán)極點,進而導致系統(tǒng)動態(tài)性能偏離預期,定量分析閉環(huán)極點隨系統(tǒng)參數(shù)變化的程度和趨勢對評價系統(tǒng)魯棒性具有較強的指導意義。

在對一般矩陣特征根求導理論[15]基礎上,本節(jié)提出一種特征根有界攝動分析方法,通過求取閉環(huán)系統(tǒng)的橫側向模態(tài)極點關于各參數(shù)的靈敏度,來定量分析閉環(huán)系統(tǒng)對參數(shù)攝動的魯棒性,以判定工程實際中控制方法可行性。

定義1 設系統(tǒng)矩陣A∈Rn×n有n個互異特征根λk(k=1,2,…,n),則稱Λk= ?λk/?A為第k個特征根λk所對應的特征根靈敏度矩陣,表示為

式中:aij為A中第i行j列的元素;Vk和Uk分別為λk所對應的右特征向量和左特征向量。

可見,特征根靈敏度矩陣表示的是矩陣中某元素對給定特征根的影響程度,該值越大,說明系統(tǒng)特征根對該參數(shù)越敏感。

式(25)推導過程如下:首先,根據(jù)特征向量定義,Vk和Uk分別滿足:

特征根靈敏度矩陣Λk體現(xiàn)了系統(tǒng)對各參數(shù)單位變化時的敏感程度,但其并不與各參數(shù)的攝動程度直接關聯(lián),為了定量體現(xiàn)不確定性對特征根的影響,本文提出了特征根有界攝動矩陣的概念。

定義2 設系統(tǒng)陣為A∈Rn×n的參數(shù)攝動量化矩陣為R,滿足:

特征根有界攝動系數(shù)為

定義2 引入了參數(shù)攝動修正項R,可以看出Γk中的每一個元素代表對應參數(shù)在攝動范圍內(nèi)對極點的最大影響,而ηk則表征所有參數(shù)最大影響的疊加,ηk越大,說明系統(tǒng)特征根的魯棒性越差,反之則魯棒性越好。

通過引入上述概念,考慮了工程應用中的實際參數(shù)攝動情況,進一步提升了其適用性,也能更準確地幫助判斷所設計的閉環(huán)系統(tǒng)對于參數(shù)攝動的魯棒性。

基于前述分析與方法,下文將結合對象模型在具體工作點給出2 種不同的控制方案,并用所提特征根有界攝動矩陣來判定2 種控制方案的魯棒性。

4 仿真結果與分析

為驗證本文方法,選取工作點為馬赫數(shù)10、高度36 km、迎角4°為例進行2 種方案的設計與分析。 首先,給出此工作點的線性化模型為

設舵偏限幅為±35°,舵機模型選擇為

可知,原系統(tǒng)荷蘭滾模態(tài)與滾轉阻尼特性都需要改善。

基于上述設定,基于六自由度非線性模型作為被控對象,本節(jié)將對所提出的2 種控制器設計方案及特征根有界攝動矩陣的計算過程進行比較和分析。

4.1 極點配置反饋控制方案

經(jīng)2.3 節(jié)分析,已經(jīng)確定為側滑角與滾轉角速率,那么主要問題在于確定反饋增益。 設期望荷蘭滾阻尼比和滾轉阻尼極點要求為

從系統(tǒng)極點來看,經(jīng)過對側滑角的大增益反饋,顯著提升了荷蘭滾阻尼,并且加大了滾轉阻尼,標稱狀態(tài)下滿足預期設計要求。 采用特征根有界攝動分析方法判斷方案1 的魯棒性,按式(33)分別求λd1、λr所對應的特征根有界攝動矩陣分別為

特征根有界攝動系數(shù)分別為

從特征根有界攝動矩陣來看,閉環(huán)系統(tǒng)極點在R所描述的攝動情況下,對L′p、N′p、L′β項的攝動比較敏感,在實際復雜環(huán)境下飛行時,可能會嚴重影響閉環(huán)極點,進而導致該控制方案魯棒性較差。

為了驗證上述結論的正確性,參照參數(shù)攝動量化矩陣(36)對L′p項加入±25%的參數(shù)攝動,即在區(qū)間[ -9.033 4, -5.420 0]變化,繪制閉環(huán)系統(tǒng)極點隨L′p變化的情況,如圖7 所示,其中線條越粗大代表L′p越大。

從圖7 中可以看出,荷蘭滾模態(tài)在加入?yún)?shù)擾動后依舊在左半平面變化,但滾轉阻尼極點變化非常劇烈,在有限的參數(shù)攝動下甚至會導致系統(tǒng)的主導極點改變,嚴重影響閉環(huán)系統(tǒng)性能。 這樣的結果與上文計算的Γr、Γd1保持了一致。

圖7 方案1 極點受L′p擾動的影響Fig.7 Poles of method 1 vary with L′p

為了進一步說明本文方法的有效性,進行了閉環(huán)六自由度非線性仿真,分別為標稱狀態(tài)和拉偏狀態(tài)仿真,參數(shù)拉偏包含氣動系數(shù)誤差、大氣密度誤差、總體參數(shù)誤差及風干擾,均取工程極端偏差值,總擾動與式(36)相匹配。 傾側角跟蹤曲線及側滑角和舵偏的響應曲線如圖8 ~圖10 所示。

圖8 方案1 傾側角響應Fig.8 Bank angle response of method 1

圖9 方案1 側滑角響應Fig.9 Sideslip angle response of method 1

圖10 方案1 舵面響應Fig.10 Aileron response of method

從仿真結果中可以看出,基于極點配置的方案在標稱狀態(tài)下,傾側角能夠穩(wěn)定跟蹤期望指令,側滑角振蕩收斂較快,系統(tǒng)舵偏響應調(diào)整迅速,除瞬態(tài)階段無明顯波動;但在參數(shù)拉偏后,控制器雖然穩(wěn)定,但傾側角出現(xiàn)較大的波動和靜差,側滑角和舵面的波動加劇且收斂變慢,系統(tǒng)性能魯棒性較差,與特征根有界攝動分析結果吻合。

4.2 反饋解耦控制方案

工作點與前文保持一致,根據(jù)式(24)反饋增益為

沿用方案1 所用R陣并按式(32)、式(33)分別求λ′d1、λ′r所對應的特征根有界攝動矩陣分別為

其對應特征根有界攝動系數(shù)為

與方案1 有顯著不同,方案2 系統(tǒng)極點總的來看對參數(shù)攝動不敏感,這說明參數(shù)擾動對系統(tǒng)極點的影響很小,系統(tǒng)魯棒性較強。 其特征根有界攝動系數(shù)要顯著小于方案1 的矩陣范數(shù),這一結果同樣印證其對參數(shù)攝動不敏感。

同樣在方案中對L′p項加入±25%的攝動,即在區(qū)間[ -7.535 0, -4.521 0]變化,系統(tǒng)極點隨L′p變化的情況如圖11 所示,其中線條越粗大代表L′p越大。

從圖11 中可以看出,采用方案2 閉環(huán)后的系統(tǒng),在加入±25% 的L′p參數(shù)攝動后,系統(tǒng)荷蘭滾極點幾乎不受影響,滾轉阻尼極點受的影響有限,魯棒性更強,這也與計算得到的Γ′r、?!鋎1結果一致。在與方案1 相同條件下進行閉環(huán)六自由度仿真,得到響應如圖12 ~圖14 所示。

圖11 方案2 極點受擾動的影響Fig.11 Poles of method 2 vary with uncertainties

圖12 方案2 傾側角響應Fig.12 Bank angle response of method 2

圖13 方案2 側滑角響應Fig.13 Sideslip angle response of method 2

圖14 方案2 舵面響應Fig.14 Aileron response of method 2

從仿真結果來看,由于偏航通道利用其本身靜穩(wěn)定性實現(xiàn)穩(wěn)定,側滑角呈現(xiàn)緩慢收斂的特點。但傾側角的控制效果良好,調(diào)節(jié)時間在2 ~3 s內(nèi),且在極端干擾情況下仍能對飛行器實現(xiàn)有效控制,相比方案1 有更強的魯棒性。

綜合來看,對極點配置方案雖然實現(xiàn)了荷蘭滾模態(tài)的快速鎮(zhèn)定,但由于其高增益特性,導致魯棒性下降,適用于模型精確情況。 而模態(tài)解耦方案的滾轉與俯仰通道控制效果更好,其優(yōu)勢在于魯棒性強,工程應用的潛力更大。

5 結 論

針對高超聲速飛行器的雙通道控制問題,本文提出了一種基于特征根有界攝動的魯棒性分析方法,并給出了2 種實現(xiàn)雙通道反饋控制的方案。

1) 給出了一種更為詳細的荷蘭滾阻尼預測式,特別是給出了交叉阻尼系數(shù)對荷蘭滾阻尼的影響,進而引出了2 種雙通道控制方案。 另外在精度上也得到了提升,全包線荷蘭滾阻尼預測精度平均提升了15%,荷蘭滾頻率預測精度平均提升了25%。

2) 所提的反饋控制魯棒性分析方法易于工程應用,引入的特征根靈敏度矩陣與特征根有界攝動矩陣計算簡單,概念本身從矩陣特征根出發(fā)也便于使用經(jīng)典控制理論研究系統(tǒng)規(guī)律。

3) 仿真表明,所提2 種控制方案可以利用副翼有效改善橫側向控制特性,實現(xiàn)雙通道控制,參數(shù)拉偏的魯棒性驗證結果與特征根有界攝動分析結果吻合。

猜你喜歡
模態(tài)
基于BERT-VGG16的多模態(tài)情感分析模型
跨模態(tài)通信理論及關鍵技術初探
一種新的基于模態(tài)信息的梁結構損傷識別方法
工程與建設(2019年1期)2019-09-03 01:12:12
多跨彈性支撐Timoshenko梁的模態(tài)分析
車輛CAE分析中自由模態(tài)和約束模態(tài)的應用與對比
國內(nèi)多模態(tài)教學研究回顧與展望
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
由單個模態(tài)構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
利用源強聲輻射模態(tài)識別噪聲源
日版《午夜兇鈴》多模態(tài)隱喻的認知研究
電影新作(2014年1期)2014-02-27 09:07:36
主站蜘蛛池模板: 午夜a级毛片| 91成人试看福利体验区| 国产美女叼嘿视频免费看| 久久中文字幕2021精品| 亚洲综合专区| 无码内射中文字幕岛国片| 国产裸舞福利在线视频合集| 毛片免费网址| 又爽又黄又无遮挡网站| 国产精彩视频在线观看| 亚洲免费人成影院| 欧美成人精品高清在线下载| 色婷婷国产精品视频| 18禁色诱爆乳网站| 永久免费AⅤ无码网站在线观看| 在线观看国产网址你懂的| 色九九视频| 视频二区欧美| 亚洲高清资源| 亚洲国产精品久久久久秋霞影院| 国产成人综合网在线观看| 老色鬼久久亚洲AV综合| 国产精品区视频中文字幕| 伊人久久青草青青综合| 欧美yw精品日本国产精品| 久久综合一个色综合网| 色综合国产| 在线无码九区| 国产成人永久免费视频| 亚洲欧洲日韩久久狠狠爱| 538国产在线| 久久五月天综合| 99久视频| 日韩免费毛片| v天堂中文在线| 伊人激情综合| 国产精品入口麻豆| 中文字幕亚洲专区第19页| 亚洲天天更新| 久久福利网| 青青草原国产一区二区| 91在线播放国产| 久久亚洲黄色视频| 亚洲欧美日韩中文字幕一区二区三区| 无码福利日韩神码福利片| 日本欧美精品| 国外欧美一区另类中文字幕| 久久婷婷色综合老司机| 浮力影院国产第一页| 免费人成在线观看成人片| 日本五区在线不卡精品| 亚洲一级毛片在线观播放| 成人午夜亚洲影视在线观看| 国产一在线| 国产高清又黄又嫩的免费视频网站| 久久国产高清视频| 国产超薄肉色丝袜网站| 爱色欧美亚洲综合图区| 国产精品无码作爱| 国产福利一区在线| 久久久久久久97| 一本一道波多野结衣一区二区| 91综合色区亚洲熟妇p| 亚洲天堂区| 欧美亚洲综合免费精品高清在线观看 | 亚洲一区二区日韩欧美gif| 99久久无色码中文字幕| 亚洲一区二区三区香蕉| 欧美激情,国产精品| 岛国精品一区免费视频在线观看| 91久久国产热精品免费| 久久中文电影| 亚洲自偷自拍另类小说| 精品人妻系列无码专区久久| 97国产精品视频自在拍| 欧美人与牲动交a欧美精品| 午夜国产不卡在线观看视频| 小说区 亚洲 自拍 另类| 广东一级毛片| 国产情侣一区| 欧美日韩中文国产| 中文天堂在线视频|