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

混合小推力航天器軌道保持高性能滑模控制

2019-08-15 03:01:56陳弈澄齊瑞云張嘉芮王煥杰
航空學(xué)報 2019年7期
關(guān)鍵詞:模型設(shè)計

陳弈澄,齊瑞云,*,張嘉芮,王煥杰

1.南京航空航天大學(xué) 自動化學(xué)院,南京 211100 2. 南京航空航天大學(xué) 先進(jìn)飛行器導(dǎo)航、控制與健康管理工業(yè)和信息化部重點實驗室,南京 211100 3. 上海航天控制技術(shù)研究所,上海 201109 4. 上海市空間智能控制技術(shù)重點實驗室,上海 201109

自1957年蘇聯(lián)發(fā)射第1顆人造衛(wèi)星以來,人類已經(jīng)開展了從近地?zé)o人、近地載人、月球無人、月球載人任務(wù)到太陽系內(nèi)各行星的探索任務(wù)[1]。深空探測使得人們更加了解地球乃至宇宙,對科技進(jìn)步和文明發(fā)展具有深遠(yuǎn)意義。然而,傳統(tǒng)的深空探測航天器越來越受到推進(jìn)火箭和化學(xué)燃料的限制,漸漸無法勝任一些太空任務(wù)。為了使深空探測的距離更遠(yuǎn),在軌完成的空間任務(wù)更多,太陽帆作為一種新型推進(jìn)方式,通過太陽光子照射到大面積輕質(zhì)薄膜帆上的反射光壓獲得動力加速度,無需消耗燃料就可以使航天器實現(xiàn)小推力連續(xù)飛行,對深空探測非常有吸引力[2]。因此,太陽帆航天器以其比沖無限大、軌道任務(wù)豐富等優(yōu)勢,得到了各國的普遍關(guān)注。2010年,日本太陽輻射驅(qū)動的行星際風(fēng)箏航天器“伊卡洛斯(IKAROS)”在世界上首次實現(xiàn)了在軌展開和性能驗證[3],更加鼓舞了人們利用太陽帆探索太空的熱情。

針對傳統(tǒng)推進(jìn)航天器的軌道控制,張海博等[4]研究了航天器系統(tǒng)的相對軌道控制問題,在推力有限情況下設(shè)計了分布式自適應(yīng)協(xié)同控制律。吳錦杰等[5]在輸入飽和條件下,進(jìn)一步研究了航天器相對運動的姿軌一體化設(shè)計問題,提出了一種魯棒自適應(yīng)控制策略。上述工作考慮到實際工程應(yīng)用中的燃料限制,均對航天器推力進(jìn)行了一定約束。針對太陽帆航天器的軌道控制,史曉寧等[6]以水星探測為背景,提出一種行星軌道分段捕獲方法并設(shè)計相應(yīng)的聯(lián)合解析最優(yōu)控制律。Farres[7]基于其在日地拉格朗日點的軌道動力學(xué)方程,通過平衡點附近穩(wěn)定流形的計算得到滿足要求的轉(zhuǎn)移軌道以及帆面指向。但總體來說,太陽帆推進(jìn)的技術(shù)成熟度仍較低,若完全依靠其執(zhí)行深空探測任務(wù)將具有很大限制。而太陽電推進(jìn)技術(shù)發(fā)展相對成熟,與化學(xué)燃料產(chǎn)生大量氣體的推進(jìn)方式相比,依靠靜電場作用加速噴出工質(zhì)電離生成的離子流,有著大幅減少推進(jìn)劑、提供精確可靠的推力矢量等優(yōu)勢,已經(jīng)被成功應(yīng)用于多項深空探測任務(wù)[8],包括JAXA的Hayabusa任務(wù)(2003),NASA的Dawn任務(wù)(2007),ESA的GOCE任務(wù)(2009)等。

混合小推力航天器是指將太陽帆技術(shù)與太陽電技術(shù)相結(jié)合應(yīng)用于推進(jìn)系統(tǒng)的連續(xù)小推力航天器,既可以解決太陽帆不能提供指向日心方向的推進(jìn)力分量的問題,又可以節(jié)省太陽電推進(jìn)劑消耗,使航天器具有更長的任務(wù)壽命。目前,針對混合小推力航天器,學(xué)者們相繼研究了地球靜止懸浮軌道[9]、高橢圓軌道[10]以及行星際轉(zhuǎn)移軌道[11-12]等非開普勒軌道任務(wù),但關(guān)于其在日心懸浮軌道的研究文獻(xiàn)還不多。

日心懸浮軌道是一種以太陽為中心,航天器推進(jìn)力與所受引力部分相平衡形成的懸浮在黃道平面上方的圓形非開普勒軌道。一方面能夠有效避免黃道平面漂浮的大量太空塵埃的影響,另一方面地球磁場以及紅外射線等帶來的干擾相對微弱,因而為在軌道上運行的航天器提供了良好的探測條件[13]。日心懸浮軌道屬于周期性軌道,當(dāng)其運轉(zhuǎn)角速度與地球公轉(zhuǎn)角速度一致時,稱其為地球同步日心懸浮軌道,由于該軌道上的航天器與地球始終保持一定方位,因此可以應(yīng)用于空間天氣監(jiān)測、極端災(zāi)害預(yù)報以及地球極區(qū)連續(xù)觀測等任務(wù)。

物理學(xué)家Forward[14]首次提出并驗證了太陽帆日心懸浮軌道。McInnes和Colin[15]通過線性化后的動力學(xué)方程求解了不同類型懸浮軌道的穩(wěn)定性,提出需要采取一定的軌道保持控制策略來保證太陽帆在目標(biāo)軌道的穩(wěn)定運行。Gong等[16]將設(shè)計問題轉(zhuǎn)化成優(yōu)化問題,通過構(gòu)造表征被動穩(wěn)定性的優(yōu)化目標(biāo)函數(shù)得到了帆面姿態(tài)參數(shù),實現(xiàn)了太陽帆在日心懸浮軌道的被動控制。Bookless和Mcinnes[17]根據(jù)忽略高階項得到的日心懸浮軌道線性化運動模型,采用最優(yōu)控制方法設(shè)計了太陽帆姿態(tài)角控制律以實現(xiàn)軌道保持控制目標(biāo)。錢航等[18]基于太陽帆航天器在參考日心懸浮軌道附近的線性化狀態(tài)方程,采用了線性二次型調(diào)節(jié)器設(shè)計了軌道保持控制律。上述針對太陽帆航天器在日心懸浮軌道的研究,大多對動力學(xué)方程進(jìn)行線性化展開,導(dǎo)致對精確模型的依賴度較高,又忽略了建模誤差、復(fù)雜深空環(huán)境等帶來的模型不確定性的影響。同時,單一太陽帆推進(jìn)存在著收斂速度慢,控制精度低的問題。張楷田等[19]在上述研究的基礎(chǔ)上,針對太陽帆、太陽電混合推進(jìn)航天器,采用了對模型依賴度較低的自抗擾方法設(shè)計軌道保持控制器,取得了較好的效果,但參數(shù)整定較為繁瑣,且出現(xiàn)了較大的超調(diào)量,即軌道位置一度偏離較大。

針對以上不足,本文根據(jù)模型不確定性的來源將其分為內(nèi)部未建模動態(tài)和外部未知擾動,設(shè)計了一種基于改進(jìn)型條件積分滑模面和徑向基(RBF)神經(jīng)網(wǎng)絡(luò)的自適應(yīng)控制方法,對于模型不確定性具有很強的魯棒性,同時具有很好的控制性能。首先,根據(jù)圓型限制性三體問題模型,建立了日心懸浮軌道柱面坐標(biāo)系下的混合小推力航天器動力學(xué)方程;其次,設(shè)計高性能滑模控制器,再將求取的虛擬控制量轉(zhuǎn)換成實際控制量;最后,數(shù)值仿真驗證了所設(shè)計控制方法的有效性。本文的主要貢獻(xiàn)有:

1) 綜合考慮模型不確定性,設(shè)計滑模控制器并采用RBF神經(jīng)網(wǎng)絡(luò)和自適應(yīng)方法進(jìn)行在線估計,提高了系統(tǒng)魯棒性。

2) 滑模控制器中基于改進(jìn)型條件積分滑模面和雙冪次趨近律進(jìn)行設(shè)計,減小了軌道位置響應(yīng)超調(diào)量,并具有較快的收斂速度。

3) 引入虛擬控制量進(jìn)行控制律的設(shè)計,再根據(jù)推進(jìn)劑最優(yōu)的條件求逆轉(zhuǎn)換成實際控制量,避免了將動力學(xué)方程線性化引起的收斂精度不高的問題。

1 混合小推力航天器日心懸浮軌道模型

1.1 太陽帆光壓推進(jìn)模型

太陽帆光壓推進(jìn)模型[20]如圖1所示,假設(shè)帆面為理想全反射平面,則光壓推進(jìn)加速度as表達(dá)式為

(1)

建立太陽帆本體坐標(biāo)系[18]Osxsyszs如圖2所示,以太陽帆質(zhì)心Os為坐標(biāo)系原點,xs軸指向日心到太陽帆位置向量r1方向,ys軸在軌道平面內(nèi)并指向太陽帆運動方向,zs軸根據(jù)右手定則確定,即指向角速度向量方向。設(shè)三軸單位向量分別為i、j、k,則根據(jù)幾何關(guān)系,n的空間方位可以表示為

n=icosα+jsinαsinδ+ksinαcosδ

(2)

式中:δ為n在Osyszs平面內(nèi)的投影與zs軸的夾角,定義為太陽帆轉(zhuǎn)角(-π≤δ≤π)。

圖1 太陽帆光壓推進(jìn)模型[20]Fig.1 Model of the solar sail propulsion[20]

圖2 太陽帆本體坐標(biāo)系[18]Fig.2 Body-fixed coordinate system of solar sail[18]

1.2 日心懸浮軌道坐標(biāo)系

建立日心懸浮軌道坐標(biāo)系[19]Oxyz如圖3所示,以日心O為坐標(biāo)系原點,x軸指向J2000星歷的春分點,z軸指向地球公轉(zhuǎn)角速度方向,y軸根據(jù)右手定則確定。設(shè)三軸單位向量分別為x0、y0、z0,則根據(jù)幾何關(guān)系,i、j、k可以表示為

(3)

航天器以角速度ω在日心懸浮軌道上運行,設(shè)(x,y,z)為航天器位置坐標(biāo),由式(1)~式(3) 可得,太陽帆光壓推進(jìn)加速度as在日心懸浮軌道坐標(biāo)系各軸上的分量為

(4)

圖3 日心懸浮軌道坐標(biāo)系[19]Fig.3 Heliocentric displaced orbit coordinate system[19]

1.3 圓型限制性三體問題模型

在太陽系運動系統(tǒng)中,當(dāng)某一天體的質(zhì)量足夠小,運動狀態(tài)不會使其所在系統(tǒng)中另外天體的運動狀態(tài)發(fā)生變化,稱作小天體,各類航天器均屬于小天體。一個小天體在另外兩個運動狀態(tài)確定的天體作用下的運動學(xué)問題稱作限制性三體問題[21],另外兩個天體仍采用二體模型,運動軌跡是以質(zhì)心為焦點的圓錐曲線。在混合小推力航天器、太陽和地球組成的三體系統(tǒng)中,將地球繞日運動軌道近似成圓軌道,研究該系統(tǒng)中航天器的運動問題就可以采用圓型限制性三體問題模型。

根據(jù)這一模型,混合小推力航天器的軌道動力學(xué)方程為

(5)

式中:μe為地心引力常數(shù);r2為地心指向航天器位置向量;ae為太陽電推進(jìn)加速度。

考慮系統(tǒng)不確定性,由來源與性質(zhì)的不同將其分為內(nèi)部未建模動態(tài)和外部未知擾動,并采用柱面坐標(biāo)系簡化分析,可以得到混合小推力航天器在日心懸浮軌道的動力學(xué)模型為

(6)

2 軌道保持控制器設(shè)計

(7)

式中:

(8)

(9)

2.1 改進(jìn)型條件積分滑模面設(shè)計

滑模控制本質(zhì)上是一種變結(jié)構(gòu)控制,所具備的不連續(xù)控制特性可以使系統(tǒng)沿指定的狀態(tài)軌跡上下作小幅、高頻的運動,即滑動模態(tài)運動。由于滑動模態(tài)可以自行設(shè)計并獨立于對象的參數(shù)及擾動,因而處于滑動模態(tài)的系統(tǒng)有著對模型不確定性不敏感的優(yōu)點[26]。滑模控制的第1步就是進(jìn)行滑模面的設(shè)計,直接決定了滑模運動的動態(tài)品質(zhì)。為了減小初始誤差較大時積分飽和引起的超調(diào)量,Seshagiri和Khalil[27]對積分項作出了調(diào)整,在傳統(tǒng)積分滑模面的基礎(chǔ)上提出了條件積分滑模面。李鵬等[28]進(jìn)一步引入調(diào)節(jié)因子,使積分削弱程度由不可調(diào)變?yōu)榭烧{(diào),提出了改進(jìn)型條件積分滑模面。

在混合小推力航天器軌道保持控制系統(tǒng)中,定義誤差變量ei(i=1,2,3)為

(10)

式中:ρd、θd和zd分別為標(biāo)稱日心懸浮軌道在柱面坐標(biāo)系下的軌道半徑、軌道角度和懸浮高度。

設(shè)計改進(jìn)型條件積分滑模面si(i=1,2,3)為

(11)

式中:ki0>0、ki1>0均為可調(diào)參數(shù),影響誤差收斂速率;μi>0為邊界層厚度;λi≥-1為積分削弱程度調(diào)節(jié)因子;sat(si/μi)為飽和函數(shù),一般定義為

(12)

以下分3種情況進(jìn)行分析:

情況1 當(dāng)|si|≤μi時

情況2 當(dāng)si>μi時

ei-1/ki0(|si|-μi)(1+λi)

情況3 當(dāng)si<-μi時

ei+1/ki0(|si|-μi)(1+λi)>ei

2.2 RBF神經(jīng)網(wǎng)絡(luò)估計器設(shè)計

(13)

(14)

式中:X為輸入向量;Γij為高斯核函數(shù)的中心向量;cij為高斯核函數(shù)的寬度參數(shù)。

(15)

2.3 自適應(yīng)滑模控制律設(shè)計

對式(11)進(jìn)行求導(dǎo),可以得到滑模面的動態(tài)方程為

(16)

將式(7)和式(13)代入,可以進(jìn)一步得到

(17)

在假設(shè)1的前提下,外部干擾d與逼近誤差ε滿足:

(18)

式中:γ>0為采用RBF神經(jīng)網(wǎng)絡(luò)進(jìn)行補償后的不確定上界。

基于雙冪次趨近律[31],設(shè)計自適應(yīng)滑模控制律為

(19)

(20)

(21)

2.4 穩(wěn)定性證明

選取Lyapunov函數(shù)為

(22)

(23)

對式(23)左右兩側(cè)積分得到

(24)

由于V(0)有界,再根據(jù)式(24)可以得到

(25)

由此推出S∈Lp(1

根據(jù)式(22)和式(24)可以得到

(26)

為了減小抖振,通常用連續(xù)的飽和函數(shù)進(jìn)一步替代控制律(19)中的符號函數(shù)[33]。

2.5 控制量轉(zhuǎn)換

通過以上設(shè)計,可以得到虛擬控制量,而實際中需要轉(zhuǎn)換成太陽帆姿態(tài)角和太陽電推進(jìn)力進(jìn)行控制。由式(9)可知,太陽電推進(jìn)加速度可以表示為

(27)

(28)

太陽電推進(jìn)劑的實時質(zhì)量變化滿足[8]:

(29)

因此,為了減少推進(jìn)劑的消耗,求解太陽帆姿態(tài)角的表達(dá)式為

(30)

式中:α*和δ*分別為實際控制輸入中的太陽帆錐角和太陽帆轉(zhuǎn)角最優(yōu)解。再根據(jù)式(27),可以得到太陽電推進(jìn)加速度在日心懸浮軌道坐標(biāo)系各軸上的加速度分量最優(yōu)解,最后通過求解實時質(zhì)量可以得到實際控制輸入中的太陽電推進(jìn)力分量Tx、Ty和Tz。

3 數(shù)值仿真

假定混合小推力航天器初始質(zhì)量為500 kg,太陽帆光壓因子為0.8,選用NEXT型號太陽電推進(jìn)器參數(shù)[34]如表2所示,假定初始入軌時軌道半徑誤差和懸浮高度誤差均為10 000 km。根據(jù)經(jīng)驗調(diào)試以及系統(tǒng)性能在各個維度的平衡,設(shè)置控制器參數(shù)如表3所示。在MATLAB R2017b環(huán)境中進(jìn)行仿真驗證,采用fmincon函數(shù)求解式(30)。

表1 標(biāo)稱日心懸浮軌道參數(shù)(無量綱)

Table 1 Parameters of nominal heliocentric displaced orbit (dimensionless)

軌道半徑懸浮高度角速度0.70.41

表2 太陽電推進(jìn)器參數(shù)Table 2 Parameters of solar electric thruster

表3 自適應(yīng)滑模控制器參數(shù)Table 3 Parameters of adaptive sliding mode controller

圖4給出了綜合考慮內(nèi)部未建模動態(tài)和外部未知擾動情況下,混合小推力航天器基于本文控制器的軌道保持響應(yīng)曲線。由圖4(a)可知,軌道半徑誤差和懸浮高度誤差在40天后收斂至2 m以內(nèi),但由于滑模控制自身的特性存在輕微抖振現(xiàn)象;由圖4(b)可知,軌道角速度在10天后穩(wěn)定為地球公轉(zhuǎn)角速度,以上結(jié)果表明本文控制器具有很強的魯棒性。圖5給出了所設(shè)計的RBF神經(jīng)網(wǎng)絡(luò)估計器對各個方向未建模動態(tài)的估計結(jié)果,可以看出估計值曲線很好地逼近了實際值曲線,精確的估計和補償對于提高系統(tǒng)控制精度具有十分重要的作用。但是,在CPU為2.8 GHz i5-8400,操作系統(tǒng)為Win 10的Dell PC機上實現(xiàn)時,每執(zhí)行一次控制算法的平均用時為0.776 ms,約為不采用RBF神經(jīng)網(wǎng)絡(luò)估計器情況下的10倍。未來還需要進(jìn)一步減少算法用時,提高控制器實時性。

圖4 日心懸浮軌道保持響應(yīng)曲線Fig.4 Response curves of station-keeping control

為了進(jìn)一步體現(xiàn)本文控制器的高性能,將軌道位置誤差分別與文獻(xiàn)[18-19]進(jìn)行比較。文獻(xiàn)[18]針對太陽帆航天器在日心懸浮軌道的線性化動力學(xué)方程,采用線性二次型調(diào)節(jié)器進(jìn)行軌道保持控制,仿真結(jié)果對比如表4所示,可以看出本文直接針對非線性動力學(xué)方程設(shè)計的滑模控制器,在更短的收斂時間內(nèi)將控制精度提升了4個數(shù)量級,同時也表明混合推進(jìn)相比于單一太陽帆推進(jìn)的控制效率更高。文獻(xiàn)[19]針對混合小推力航天器,采用對精確模型依賴度較低的自抗擾控制方法,仿真結(jié)果如圖6所示,可以看出相對位置誤差在90天后收斂至50 m以內(nèi),但軌道位置一度偏離較大。而本文基于改進(jìn)型條件積分滑模面設(shè)計的控制器,可以有效地減小軌道位置超調(diào),并且具有很小的穩(wěn)態(tài)誤差。

圖7給出了軌道保持控制量曲線,由圖7(a)可知,存在初始入軌誤差的情況下,為了盡快運行在標(biāo)稱軌道,太陽帆姿態(tài)角在1天內(nèi)具有較大的變化,10天后趨于穩(wěn)定,錐角保持在23°,而轉(zhuǎn)角一直保持在0°,驗證了懸浮軌道對太陽帆帆面指向的要求,即法向量需要位于日心到太陽帆位置向量和運行角速度向量構(gòu)成的平面內(nèi);由圖7(b)可知,太陽電推進(jìn)力經(jīng)過10天后保持在6 mN以內(nèi),數(shù)值上的變化一方面用來提供航天器在懸浮軌道進(jìn)行圓周運動所需的向心力,另一方面用來克服時變的內(nèi)部未建模動態(tài)與外部擾動給系統(tǒng)帶來的影響。同時,圖8給出了不同推進(jìn)方式下的航天器質(zhì)量變化曲線,可以看出經(jīng)過1年時間,基于混合推進(jìn)僅消耗8.09 kg推進(jìn)劑,而基于單一太陽電推進(jìn)需要消耗77.94 kg推進(jìn)劑,前者可以節(jié)省約89.6%推進(jìn)劑,表明混合推進(jìn)在長周期任務(wù)上具備很大優(yōu)勢。

圖5 未建模動態(tài)及其估計曲線Fig.5 Curves of unmodeled dynamics and its estimation

表4 軌道位置誤差對比Table 4 Comparison of tracking position errors

控制器收斂時間/d軌道半徑誤差/km懸浮高度誤差/km本文40-0.00120.0016文獻(xiàn)[18]730-30.0692-198.9652

圖6 文獻(xiàn)[19]軌道位置誤差響應(yīng)曲線Fig.6 Response curves of tracking errors in Ref.[19]

圖7 日心懸浮軌道保持控制量曲線Fig.7 Control variable curves of station-keeping control

圖8 航天器質(zhì)量變化曲線Fig.8 Mass curves of spacecraft

4 結(jié) 論

1) 本文基于改進(jìn)型條件積分滑模面和雙冪次趨近律理論設(shè)計的高性能軌道保持控制器,有效減小了軌道位置超調(diào),軌道跟蹤速度快,控制精度高,對內(nèi)部未建模動態(tài)和外部未知擾動均具有較強的魯棒性。

2) 所設(shè)計的RBF神經(jīng)網(wǎng)絡(luò)對各個方向的未建模動態(tài)能夠進(jìn)行有效的在線估計和補償,將神經(jīng)網(wǎng)絡(luò)中的最優(yōu)權(quán)值向量作為不確定參數(shù)進(jìn)行自適應(yīng)調(diào)節(jié),保證了誤差的收斂性和估計的實時性。

3) 引入的虛擬控制量根據(jù)推進(jìn)劑最優(yōu)條件求逆轉(zhuǎn)換成實際控制量,既解決了將動力學(xué)方程線性化帶來的控制精度不高的問題,又體現(xiàn)了混合推進(jìn)兼具太陽帆推進(jìn)和太陽電推進(jìn)的優(yōu)勢,相比于單一太陽帆推進(jìn)控制效率更高,相比于單一太陽電推進(jìn)控制消耗更低,因此混合小推力航天器更適合長期執(zhí)行復(fù)雜的深空探測任務(wù)。

猜你喜歡
模型設(shè)計
一半模型
重要模型『一線三等角』
何為設(shè)計的守護(hù)之道?
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
《豐收的喜悅展示設(shè)計》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設(shè)計萌到家
設(shè)計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設(shè)計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 激情国产精品一区| 日韩中文无码av超清| 欧美午夜在线播放| 精品国产免费观看| 亚洲中文无码av永久伊人| 亚洲精品少妇熟女| 无码一区二区三区视频在线播放| 精品一区二区无码av| 一级毛片免费不卡在线| 久久久久九九精品影院 | 97无码免费人妻超级碰碰碰| 无码国产偷倩在线播放老年人| 国产成人福利在线视老湿机| 中文字幕亚洲另类天堂| 国产精品男人的天堂| 精品欧美一区二区三区在线| 国产91在线|日本| 伊人欧美在线| 欧美日韩国产综合视频在线观看| 国产精品中文免费福利| 8090午夜无码专区| 久久国语对白| 欧美人与牲动交a欧美精品| 四虎综合网| 中文字幕在线日本| 国产一区二区三区视频| 国产99免费视频| 亚洲日本中文综合在线| 亚洲女同欧美在线| 欧美一级夜夜爽www| 午夜啪啪福利| 亚洲精品国产精品乱码不卞| 熟女成人国产精品视频| 亚洲国产精品国自产拍A| 国产91色| 国产a v无码专区亚洲av| 麻豆精选在线| 国产一区二区人大臿蕉香蕉| 国产精品久久久久婷婷五月| 99久久成人国产精品免费| 粉嫩国产白浆在线观看| 亚洲精品日产AⅤ| 青青极品在线| 亚洲一级毛片免费观看| 久草性视频| 欧美日韩中文字幕二区三区| 国产成人亚洲毛片| 成年免费在线观看| 亚洲综合二区| 国产激情国语对白普通话| 国产产在线精品亚洲aavv| 日韩av无码DVD| 欧美亚洲国产精品久久蜜芽| 欧美一级高清片欧美国产欧美| 在线精品亚洲国产| 91美女在线| 婷婷色狠狠干| 国产性猛交XXXX免费看| 婷婷色狠狠干| 成人91在线| 污网站在线观看视频| av一区二区三区在线观看| 国产男女免费视频| av一区二区三区在线观看| 欧美一级片在线| 久久精品91麻豆| 波多野结衣中文字幕一区二区| 最新国产麻豆aⅴ精品无| 91精品aⅴ无码中文字字幕蜜桃| 久久综合AV免费观看| 一区二区三区四区精品视频| 国产日韩精品欧美一区喷| 亚洲第一视频网| 91久久国产热精品免费| 国产精品手机在线播放| 女人18一级毛片免费观看| 亚洲欧美成人综合| 国产综合日韩另类一区二区| 91啪在线| 国产原创演绎剧情有字幕的| 高清久久精品亚洲日韩Av| 在线精品自拍|