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

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx

高超聲速飛行試驗助推段彈道優(yōu)化方法

2016-10-20 10:53:29張子明王西耀
彈道學報 2016年3期
關鍵詞:優(yōu)化模型

張子明,蔣 勁,王西耀

(中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所 高超聲速沖壓發(fā)動機技術重點實驗室,四川 綿陽 621000)

?

高超聲速飛行試驗助推段彈道優(yōu)化方法

張子明,蔣勁,王西耀

(中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所 高超聲速沖壓發(fā)動機技術重點實驗室,四川 綿陽 621000)

為減小高超聲速飛行試驗助推段飛行環(huán)境的惡劣程度,提出了一種生成三自由度最優(yōu)彈道的方法。采用基于攻角編碼的遺傳算法,并用三次樣條對攻角-時間歷程進行平滑處理。優(yōu)化過程首先搜索滿足窗口參數(shù)的可行解,然后以攻角范圍、最大動壓、最大法向過載等參數(shù)為目標函數(shù),以窗口參數(shù)為約束,搜索飛行環(huán)境惡劣程度最低的最優(yōu)化彈道。計算結果表明,該算法能夠避免理論彈道中存在局部高風險階段的情況,并具有良好的魯棒性。

高超聲速飛行器;理論彈道;遺傳算法;多目標優(yōu)化

高超聲速技術是當今世界各國競相研究的熱點。該技術由于涉及激波/邊界層干擾、湍流轉捩、非平衡流動等許多復雜的空氣動力學關鍵問題,國內(nèi)外仍在攻關之中,普遍采用地面試驗、數(shù)值計算和飛行試驗相結合的研究手段。

模型飛行試驗[1-2]是針對單項或者多項空氣動力學關鍵技術的飛行試驗。由于試驗模型在高馬赫數(shù)下才能開始工作,因此,需要采用助推的方式,將其送入預定的試驗窗口。在助推段,重力和負升力的綜合作用使試驗模型進入水平飛行的試驗窗口。本文研究的飛行試驗模型采用類似于圖1所示的升力體外形[3]。與其他飛行試驗相比,高超聲速試驗模型的飛行試驗具有以下特點:①升力體外形的試驗模型壓心靠前、質(zhì)心靠后,導致舵面配平性能較差,只能在小攻角條件下飛行;②試驗模型的尖銳前緣對熱環(huán)境提出了更高的要求,需要盡量降低助推段的飛行動壓;③受試驗模型質(zhì)量及結構的約束,飛行過程中的法向過載不宜太大;④跨聲速段為高風險飛行階段,應避免大攻角飛行;⑤試驗窗口的彈道參數(shù)通常較為苛刻,允許的偏差較小。

圖1 飛行試驗模型示意圖

基于上述特點,高超聲速試驗模型的助推段彈道設計應考慮以下4個關鍵因素:①攻角變化范圍。針對升力體構型的試驗模型,大的負攻角狀態(tài)下縱向氣動特性可能是靜不穩(wěn)定的,當與某些側滑角組合時,航向也是靜不穩(wěn)定的。另外,強烈的縱橫耦合力矩和上、下氣動舵面效率差異顯著的問題,也將給橫航向控制帶來難度。因此,彈道設計應避免采用大的負攻角。②跨聲速段攻角限制。試驗模型的空氣動力非線性特征突出,跨聲速段氣動力和力矩特性變化劇烈,為高風險飛行階段,應當避免采用大的彈道機動。例如,美國飛馬座空射運載火箭以Ma=0.8投放后,以零攻角加速(同時自然降高)到約385 m/s(Ma=1.3)后才開始以大攻角向上拉升。③舵偏角變化范圍。模型的升降舵效率相對較低,導致升降舵偏角較大,且在拉偏情況下,由于舵偏角與舵偏效率的非線性變化規(guī)律容易導致舵偏飽和的情況,因此,需要盡量降低理論彈道的配平舵偏角。④進入窗口前的姿態(tài)調(diào)整。考慮到助推器的拖尾段不利于對攻角的精確控制和穩(wěn)定,且助推器工作時間存在一定的不確定量,因此,彈道設計時,在進入窗口前需要盡量減小攻角變化率,并提前將攻角調(diào)整到零攻角附近。

由此可見,高超聲速飛行試驗的助推段彈道優(yōu)化,需要針對由特定的助推裝置、氣動性能、操穩(wěn)特性、質(zhì)量特性及結構約束等環(huán)節(jié)構成的飛行系統(tǒng),尋找與該系統(tǒng)最匹配的飛行彈道,讓試驗模型以最“舒服”的方式飛行到試驗窗口,并盡量降低飛行環(huán)境的惡劣程度,確保試驗成功。這一優(yōu)化目標與常規(guī)的能量最小或時間最短等優(yōu)化目標存在顯著的區(qū)別。

目前,關于優(yōu)化計算的方法和理論已經(jīng)較為全面[3],這些最優(yōu)化理論提供了多種優(yōu)化的方法。常用的優(yōu)化方法有動態(tài)規(guī)劃、二次規(guī)劃法、遺傳算法等方法。

動態(tài)規(guī)劃是美國數(shù)學家Bellman R等人于1951年提出的解決多階段決策過程的最優(yōu)化方法。該方法在工程技術、經(jīng)濟、工業(yè)生產(chǎn)及軍事等部門都有廣泛的應用。但是,動態(tài)規(guī)劃方法也存在兩大弱點[4-6]:一是利用最優(yōu)原理得到函數(shù)方程后,沒有一種統(tǒng)一的處理方法,需要根據(jù)具體問題結合數(shù)學技巧來求解;二是所謂的“維數(shù)障礙”,即當變量個數(shù)太大時,受計算機限制而無法求解。例如,Harada A,Miyazawa Y等人采用動態(tài)規(guī)劃法對民航客機的縱向軌道進行了優(yōu)化[7-8]。盡管這些研究工作已經(jīng)極力考慮了“維數(shù)障礙”的問題,但當變量數(shù)達到20以上時,計算依然顯得困難。

二次規(guī)劃已經(jīng)成為運籌學、經(jīng)濟數(shù)學、系統(tǒng)分析和組合優(yōu)化科學的基本方法。在航空航天領域,Ran D采用序列二次規(guī)劃算法,研究了最小時間攔截問題[9],鄭總準等人采用序列二次規(guī)劃算法,對再入軌道進行了優(yōu)化[10],其最優(yōu)化目標為總加熱量最小。劉欣等人采用序列二次規(guī)劃法對助推導彈進行了彈道優(yōu)化[11],傅瑜等人也采用序列二次規(guī)劃法,對可重復使用助推器飛行器的軌跡進行了優(yōu)化研究[12]。

遺傳算法(GA)是Holland J于1975年提出的一種智能的全局優(yōu)化算法[13]。近年來,采用遺傳算法進行彈道優(yōu)化的研究也不少。例如,齊照輝等人采用遺傳算法對彈道導彈助推段軌道進行了優(yōu)化[14],楊遵等人采用遺傳算子粒子群優(yōu)化算法進行任務優(yōu)化分配,對無人機的航路進行了規(guī)劃[15]。

所有上述的彈道優(yōu)化工作,都是值得肯定和借鑒的。但是,高超聲速飛行試驗助推段的彈道優(yōu)化,需要更多考慮的是在滿足試驗窗口的前提下盡量降低飛行環(huán)境的惡劣程度,盡量降低飛行過程中的最大動壓、最大法向過載、攻角范圍等參數(shù),而不是單純的能量最優(yōu)、時間最短或總加熱量最小等單一的優(yōu)化目標。這樣的優(yōu)化目標在以前的彈道設計中不曾提及,因此,需要研究這種特殊的彈道優(yōu)化問題。

1 彈道優(yōu)化方法

1.1運動方程及優(yōu)化目標

本文采用鉛垂面內(nèi)的三自由度模型,運動方程組為[16-17]

(1)

(2)

(3)

式中:FP為助推器推力;FX,FY分別為氣動阻力和升力;α為飛行攻角;θ為彈道傾角;MZ為氣動質(zhì)心矩;ωZ為繞質(zhì)心的俯仰角速度。

終端條件為助推段結束時分離時刻的試驗窗口參數(shù):

hend=hq+hΔ

(4)

Maend=Maq+MaΔ

(5)

θend=θq+θΔ

(6)

αend=0+αΔ

(7)

(8)

以上為理論彈道的5個終端條件,其中,h為飛行高度,Ma為飛行馬赫數(shù),下標“Δ”表示允許的偏差,下標“end”表示助推段結束時刻的彈道參數(shù),下標“q”表示試驗窗口要求的期望值。其中,試驗模型質(zhì)量參數(shù)、試驗模型氣動性能、助推器推力特性為已知條件。優(yōu)化目標可以描述為

ny,max|0≤t≤tend=ny,min

(9)

pmax|0≤t≤tend=pmin

(10)

|α|max|0≤t≤tend=αmin

(11)

(12)

|α|max|0.8≤Ma≤1.5=αk,min

(13)

1.2優(yōu)化算法

本文采用遺傳算法對上述三自由度彈道進行優(yōu)化。

1.2.1編碼

編碼是遺傳算法的基礎,本文提出了一種有理數(shù)編碼方案。在助推段上按順序選取若干時間節(jié)點并組成一維數(shù)組t=(t1,t2,…,tn),其中t1=0,tn=tend,選取由每個節(jié)點時刻的飛行攻角值構成的一維數(shù)組α=(α1,α2,…,αn)為彈道個體的染色體,即基因序列。采用這種編碼方案的理由是:

①終端條件方程(7)可以自動滿足;

②可以直接判斷一部分針對攻角的約束條件,若不滿足約束,不必計算該條彈道,從而很簡單地避免了“攻角超界”的情況;

③在三自由度彈道中,攻角為物理意義上的自變量,當攻角歷程確定后可以直接得到完整的彈道,若選取其他參數(shù)作為控制對象,比如彈道傾角或飛行高度等參數(shù),都需要通過迭代計算出攻角值,無疑增加了計算復雜度。

④在工程應用中很少選用攻角作為控制對象,其主要原因在于干擾的存在和攻角測量的精度不足。但理論彈道設計不存在這樣的情況。

1.2.2插值函數(shù)

僅有節(jié)點時刻的攻角值是不夠的,需要進行插值計算。插值的方法很多,通常有線性插值、多項式插值、樣條插值等方法。本文采用三次樣條插值,理由如下:

①終端條件方程(8)可以自動滿足;

②算法易于實現(xiàn),當節(jié)點個數(shù)較多時,采用多項式插值時,方程的求解較為困難;

③通過樣條插值得到的攻角歷程是一條光滑的曲線,而不是線性插值得到的一條折線,這更有利于后期的飛行控制系統(tǒng)設計工作。

1.2.3適應值函數(shù)

適應值是對染色體(個體)進行評價的一種指標,是GA進行優(yōu)化作用的主要信息。在生物界中,對生存環(huán)境適應程度較高的物種或個體將有更多的繁殖機會,反之則機會較少,甚至滅絕。與之類似,遺傳算法也使用適應度這個概念來量度種群中各個個體在優(yōu)化計算中有可能達到或接近于或有助于找到最優(yōu)解的優(yōu)良程度。

適應值函數(shù)通常是由優(yōu)化的目標函數(shù)變換而成的。本文提出了一種將優(yōu)化過程劃分為2個階段的方法:首先搜索滿足終端條件的可行解,然后在可行解的基礎上搜索飛行環(huán)境惡劣程度更小的彈道。2個階段的適應值函數(shù)分別為

A1=|hend-hq|k1+|Maend-Maq|k2+|θend-θq|k3

(14)

(15)

式中:k1,k2,…,k8為人為設定的正系數(shù)。針對第2階段的適應值函數(shù)A2,需要增加懲罰函數(shù),確保滿足試驗窗口條件。值得一提的是,上述適應值函數(shù)應取極小值為優(yōu)化目標。

在實際計算過程中,系數(shù)的確定還需要考慮各項參數(shù)本身的物理意義。比如,若分離點高度與終端條件偏差10 m,或者動壓偏差1 kPa,是完全可以接受的,但是彈道傾角與窗口條件即使僅偏差0.1 rad,就無法接受了。在確定以上系數(shù)時應充分考慮這一點。

2 算例

本文研究的算例采用固體火箭發(fā)動機為助推動力,從試驗場地垂直發(fā)射(海拔高度約為1 km),并要求助推結束時進入飛行試驗窗口。助推段飛行時間約為86 s,選取時間節(jié)點的間隔為2 s,即扣除初始和終端時刻的攻角為已知后,實際參與編碼的基因為42個。采用的種群大小為100,選擇概率為0.9,并保留5%的最優(yōu)個體,變異概率為0.05,種群初值為全程零攻角飛行。圖2~圖7給出了最終彈道參數(shù)曲線。

圖2 攻角隨時間變化曲線

從圖2可知,試驗模型的助推段全程為負攻角飛行。在助推段后期,助推火箭發(fā)動機進入拖尾段之前的一段較長時間范圍內(nèi),試驗模型保持-5.3°左右的最大負攻角飛行。這樣的飛行方式既充分利用負升力來調(diào)整飛行的彈道傾角,同時又避免了局部出現(xiàn)大攻角的情況,從而降低了最大法向過載和縱向配平舵偏,達到了減小飛行環(huán)境惡劣程度的優(yōu)化目標。助推器熄火之前8 s,開始將攻角逐漸調(diào)整為0°以滿足頭體分離對姿態(tài)的要求。從圖3可知,跨聲速段出現(xiàn)在彈道的第23 s到第36 s之間,期間的攻角約為-1.6°。圖4給出了飛行高度h及縱向配平舵偏角δZ隨時間的變化曲線。由圖可知:助推結束時,飛行高度約為25 km;最大舵偏角約為4°。后期的仿真分析表明,全助推段彈道可控,且有較大的控制裕度。

圖3 飛行馬赫數(shù)隨時間變化曲線

圖4 飛行高度及縱向配平舵偏角隨時間變化曲線

圖5為動壓p和法向過載ny隨時間的變化規(guī)律。由圖可知,最大動壓約為100 kPa,且在助推段出現(xiàn)了2次,分別是低馬赫數(shù)低空階段和高馬赫數(shù)高空階段。在通常情況下,助推彈道動壓曲線的這2個峰值是不相等的。其中,由于受到試驗窗口參數(shù)的限制,高馬赫數(shù)高空階段的動壓峰值變化范圍不會太大。當?shù)婉R赫數(shù)低空階段的動壓峰值大于高馬赫數(shù)高空階段的動壓峰值時,需要調(diào)整助推段的攻角歷程,以減小低馬赫數(shù)低空階段的動壓峰值。調(diào)整的結果導致彈道在40 s附近時出現(xiàn)攻角曲線的微小波動。這種現(xiàn)象避免了局部出現(xiàn)高動壓的情況,體現(xiàn)了彈道優(yōu)化的結果。法向最大過載約為2.2g,出現(xiàn)在進入試驗窗口前即將開始調(diào)整姿態(tài)的時刻。

圖5 飛行動壓及法向過載隨時間變化曲線

圖6給出了改變種群的初始值對計算結果的影響。將各個節(jié)點的初始攻角統(tǒng)一設為±1°以及±2°,該優(yōu)化算法均能最終達到與上述彈道基本一致的結果,這表明本文采用的優(yōu)化算法并不依賴于計算初值,算法具有較好的魯棒性。

圖6 不同初值的攻角計算結果

3 結論

本文采用遺傳算法對高超聲速飛行試驗的助推段彈道進行了優(yōu)化。結論如下:

①以試驗窗口參數(shù)為約束條件,以降低飛行環(huán)境的惡劣程度為彈道優(yōu)化的目標,可以既滿足試驗窗口條件,又最大限度地降低助推段的最大攻角、最大飛行動壓、最大法向過載,從而降低氣動布局、操穩(wěn)特性、結構強度等分系統(tǒng)的設計難度。

②采用該算法得到的理論彈道,結果合理,且算法不依賴計算初值,具有良好的魯棒性。

[1]甘楚雄,劉冀湘.彈道導彈與運載火箭總體設計[M].北京:國防工業(yè)出版社,1996.

GAN Chu-xiong,LIU Ji-xiang.Ballistic missiles and launch vehicles design[M].Beijing:National Defense Industry Press,1996.(in Chinese)

[2]龍樂豪.導彈與航天叢書:總體設計[M].北京:中國宇航出版社,1991.

LONG Le-hao.Missile and space series:overall design[M].Beijing:China Astronautic Publishing House,1991.(in Chinese)

[3]JOSEPH M H.The X-51A scramjet engine flight demonstration program,AIAA 2008-2 540[R].2008.

[4]何獻忠,李萍.優(yōu)化技術及其應用[M].北京:北京理工大學出版社,1995:12-156.

HE Xian-zhong,LI Ping.Optimization technique and applications[M].Beijing:Beijing Institute of Technology Press,1995:12-156.(in Chinese)

[5]方洋旺.隨機系統(tǒng)最優(yōu)控制[M].北京:清華大學出版社,2005.

FANG Yang-wang.Stochastic system optimal control[M].Beijing:Tsinghua University Press,2005.(in Chinese)

[6]張光澄.非線性最優(yōu)化計算方法[M].北京:高等教育出版社,2003.ZHANG Guang-cheng.Computational methods for nonlinear optimization[M].Beijing:High Education Press,2003.(in Chinese)

[7]HARADA A.Dynamic programming trajectory optimization by piecewise linear approximation,AIAA 2015-1 075[R].2015.

[8]MIYAZAWA Y.Dynamic programming application to airliner four dimensional optimal flight trajectory,AIAA 2013-4 969[R].2013.

[9]DAI R.Three-dimensional minimum-time interception trajectory planning using nonlinear programming and collocation,AIAA 2007-6 787[R].2007.

[10]鄭總準,吳浩,王永驥.基于序列二次規(guī)劃算法的再入軌跡優(yōu)化研究[J].航天控制,2009,27(6):8-13.

ZHENG Zong-zhun,WU Hao,WANG Yong-ji.Reentry trajectory optimization using sequential quadratic programming[J].Aerospace Control,2009,27(6):8-13.(in Chinese)

[11]劉欣,楊濤,張青斌.助推-滑翔導彈彈道優(yōu)化與總體參數(shù)分析[J].彈道學報,2012,24(1):43-48.

LIU Xin,YANG Tao,ZHANG Qing-bin.Trajectory optimization and parameter analysis for boost-glide missile[J].Journal of Ballistics,2012,24(1):43-48.(in Chinese)

[12]傅瑜,馬永青,崔乃剛.可重復使用助推飛行器軌跡優(yōu)化研究[J].彈道學報,2011,23(4):5-9.

FU Yu,MA Yong-qing,CUI Nai-gang.Research on optimal trajectory of reusable boost vehicle[J].Journal of Ballistics,2011,23(4):5-9(in Chinese)

[13]王凌.智能優(yōu)化算法及其應用[M].北京:清華大學出版社,2001:36-59.WANG Ling.Intelligent optimization algorithms with applications[M].Beijing:Tsinghua University Press,2001:36-59.(in Chinese)

[14]齊照輝,梁偉,劉雪,等.彈道導彈助推段突防多目標優(yōu)化設計方法[J].彈道學報,2010,22(2):39-43.

QI Zhao-hui,LIANG Wei,LIU Xue,et al.Multiobjective optimization algorithm of survivability of ballistic missile in boost phase[J].Journal of Ballistics,2010,22(2):39-43.(in Chinese)

[15]楊遵,雷虎民,唐治理.基于作戰(zhàn)效能的無人機航路規(guī)劃研究[J].電光與控制,2009,16(11):14-18.

YANG Zun,LEI Hu-min,TANG ZHI-li.Combat effectiveness based path planning for unmanned air vehicles[J].Electron Optics & Control,2009,16(11):14-18.(in Chinese)

[16]韓子鵬.彈箭外彈道學[M].北京:北京理工大學出版社,2008.

HAN Zi-peng.Exterior ballistics of projectile and rocket[M].Beijing:Beijing Institute of Technology Press,2008.(in Chinese)

[17]王中原,周衛(wèi)平.外彈道設計理論與方法[M].北京:科學出版社,2004:59-102.

WANG Zhong-yuan,ZHOU Wei-ping.Theory and method of exterior ballistic design[M].Beijing:Science Press,2004:59-102.(in Chinese)

Investigation on Boost Trajectory Optimization Method of Hypersonic Flight Test

ZHANG Zi-ming,JIANG Jin,WANG Xi-yao

(Science and Technology on Scramjet Laboratory,Hypervelocity Aerodynamics Institute,CARDC,Mianyang 621000,China)

To reduce the badness of the flight environment of the boost phase in the hypersonic flight test,a novel optimization method to generate 3-DOF optimal trajectory was proposed.The method based on the angle of attack(AOA)encoding genetic algorithm was adopted,and the cubic spline was employed to smooth attack angle time courses.During the optimization process,the feasible solution was first searched to meet window parameters.Then the optimal trajectory with the lowest badness of flight environment was searched for the objective functions such as attack angle range,the maximum dynamic pressure and the maximum normal overload parameters under the constraint of the windows parameters.The results show that the algorithm can avoid the local high risk in the theoretical trajectory,and it has good robustness.

hypersonic vehicle;theoretical trajectory;genetic algorithm;multi-objective optimization

2016-03-05

張子明(1974- ),男,高級工程師,研究方向為飛行試驗技術。E-mail:fatmanzzm@sina.com。

TJ013

A

1004-499X(2016)03-0012-05

猜你喜歡
優(yōu)化模型
一半模型
超限高層建筑結構設計與優(yōu)化思考
民用建筑防煙排煙設計優(yōu)化探討
關于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
由“形”啟“數(shù)”優(yōu)化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
404 Not Found

404 Not Found


nginx
主站蜘蛛池模板: 97精品久久久大香线焦| 午夜福利网址| 国产产在线精品亚洲aavv| 国产91特黄特色A级毛片| 欧美有码在线| 老司国产精品视频| 亚洲天堂视频在线播放| 欧美亚洲一二三区| 亚洲精品成人福利在线电影| 国产高清在线观看91精品| 亚洲啪啪网| 国产又黄又硬又粗| 国产jizzjizz视频| 日本91在线| 一本一道波多野结衣一区二区 | 国产精彩视频在线观看| 乱系列中文字幕在线视频 | 欧美亚洲日韩中文| 四虎国产永久在线观看| 国产主播一区二区三区| 国产系列在线| 一级毛片免费不卡在线视频| 久久久久88色偷偷| 夜夜爽免费视频| 日韩毛片基地| 国产激情无码一区二区免费| 国产乱人伦精品一区二区| 成色7777精品在线| 成人免费一级片| 一区二区三区国产| 免费精品一区二区h| 成人福利在线视频免费观看| 91久久天天躁狠狠躁夜夜| 色精品视频| 久久人人97超碰人人澡爱香蕉| 亚洲国产综合精品中文第一| 99爱视频精品免视看| 日本伊人色综合网| 青青草原国产免费av观看| 视频二区亚洲精品| 亚洲日韩第九十九页| 国内黄色精品| 国产91小视频| 伊人天堂网| 亚洲高清免费在线观看| 亚洲欧美不卡视频| 日韩精品欧美国产在线| 国产成人综合亚洲欧美在| 精品视频福利| 国产一区二区免费播放| 小13箩利洗澡无码视频免费网站| 午夜a视频| 久久久亚洲色| 亚洲高清中文字幕在线看不卡| 四虎国产在线观看| 国产精品妖精视频| 九九久久99精品| 极品性荡少妇一区二区色欲| 911亚洲精品| 欧美一级黄片一区2区| 日韩国产综合精选| 欧美日韩资源| 88av在线| 日本草草视频在线观看| 黄色网页在线播放| 亚洲精品国产综合99久久夜夜嗨| 国产高清毛片| 毛片网站在线看| 日韩经典精品无码一区二区| 欧美一级高清片欧美国产欧美| 日韩AV无码一区| 在线免费a视频| 成人国产精品一级毛片天堂| 国产精品部在线观看| 国产欧美一区二区三区视频在线观看| 国产激情第一页| 亚洲男人的天堂在线| 99精品热视频这里只有精品7| 精品国产美女福到在线不卡f| 免费观看精品视频999| 色偷偷一区二区三区| 一本大道无码高清|