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

載人登月著陸器奔月窗口搜索方法

2017-11-17 10:23:31賀波勇李海陽
航空學報 2017年4期
關鍵詞:模型

賀波勇, 李海陽

國防科學技術大學 航天科學與工程學院, 長沙 410073

載人登月著陸器奔月窗口搜索方法

賀波勇, 李海陽*

國防科學技術大學 航天科學與工程學院, 長沙 410073

對環月軌道共面交會的載人登月任務中,著陸器(LM)奔月零窗口與軌道參數精確快速設計方法進行了研究。任務采用人貨分離奔月模式,著陸器于載人飛船到達環月軌道前抵達環月共面交會軌道,著陸器近月點一次共面減速完成近月制動。提出一種三層快速精確奔月窗口搜索方法:第一層采用地心二體軌道理論解析計算月窗口及奔月軌道參數初值,作為正確性基本參考;第二層采用改進的雙二體解析動力學模型求解月窗口內奔月軌道參數變化規律;第三層采用高精度軌道動力學模型和SQP_Snopt優化求解奔月零窗口及軌道參數精確解。仿真結果表明,本文提出的三層逐級奔月窗口搜索方法能快速精確求解載人登月任務中著陸器奔月窗口及精確軌道參數,也揭示了影響著陸器奔月窗口的主次因素和規律,為中國未來載人登月工程提供參考。

載人登月; 著陸器(LM); 環月交會; 奔月窗口; 逐級搜索

月球是地球唯一的天然衛星,是人類探索地球以遠的必經之路。美國2009年6月18日發射的“月球偵查軌道器(Lunar Reconnaissance Orbiter,LRO)/月坑觀測與探測衛星(Lunar CRater Observation and Sensing Satellite,LCROSS)”,及其隨行的半人馬座火箭上面級,成功對月球南極進行了兩次撞擊,發現水冰存在的證據[1-2],引發了生命科學、天體物理等領域科學家關注,又一次將載人月球探測升溫[3]。

受重型火箭運載能力約束及現代載人航天高可靠性要求,采用人貨分離奔月及環月軌道交會的載人登月飛行模式[4]隨著中國交會對接技術成熟應運而生[5]。該飛行模式中,載人飛船通常采用地月轉移約3 d的自由返回軌道或混合軌道,而著陸器不受航天員生保系統時間約束,一般采用地月轉移能量較低的5 d飛行軌道,要求近地單脈沖共面加速出發,近月點單脈沖共面制動,之后只進行近月點高度維持,等待環月共面交會。鄭愛武和周建平[6-7]概述了載人登月自由返回軌道、混合軌道及任務中止應急返回軌道的一般設計思路及約束條件,并基于Lambert算法快速設計月地轉移軌道窗口,得到相關窗口與最小再入角基本關系。周文艷和楊維廉[8]在“嫦娥工程”中研究了月球探測器地月轉移軌道設計方法。Peng和Shen等[9-10]將雙二體假設模型解作為高精度動力學模型迭代求解初值,采用粒子群優化(Particle Swarm Optimization,PSO)算法和Multi-Start收斂算法,分別求解了繞月自由返回軌道特性和定點返回軌道特性。Li等[11]采用雙二體假設仿真求解了月地轉移窗口及軌道參數特性。賀波勇等[12]采用改進微分校正方法快速求解地月轉移軌道,但未給出結合工程實際約束的窗口求解方法。張洪禮和羅欽欽等[13-14]嘗試采用無跡卡爾曼濾波(Unscented Kalman Filter,UKF)參數估計的三體Lambert算法求解地月引力輔助變軌問題,雖然已有技術突破,但離工程使用尚有一定距離。綜上所述,目前關于月球探測器奔月任務的研究局限于兩個方向:一是單純的軌道設計方法研究;二是單純的基于雙二體假設的月窗口區間求解方法研究。缺少結合具體復雜任務(如載人登月)約束條件的地月轉移零窗口與精確軌道耦合設計方法。

本文主要研究環月軌道共面交會的載人登月任務中,著陸器奔月零窗口與軌道參數精確快速設計問題。提出了三層逐級求解策略:第一層采用地心二體模型解析估計月窗口范圍;第二層采用改進的雙二體假設解析模型快速打靶求解月窗口范圍內參數變化規律;第三層采用SQP_Snopt算法[15]和高精度軌道動力學模型求解精確零窗口及軌道參數。數值仿真算例驗證了本文方法的快速精確性,并給出了商業軟件STK的驗證結果。

1 問題描述與動力學模型

1.1 問題描述

人貨分離奔月及環月軌道交會的載人登月飛行模式可用如圖 1 所示甘特圖描述。圖中:LEO(Low Earth Orbit)為近地停泊軌道;LLO(Low Lunar Orbit)為環月圓軌道;LM(Lunar Module)為著陸器;CSM(Command and Service Module)為指揮服務推進飛行器;CEV(Crew Exploration Vehicle)為乘員探索飛行器(載人飛船);RVD(RendezVous and Docking)為交會對接;飛行器A代表著陸器,分為A1下降級和A2上升級;B為載人飛船,分為B1服務艙和B2返回艙;時刻均用T表示,上標A/B區別飛行器,下標0,1,2,…代表任務主要窗口或軌道機動時刻序列。

圖1 人貨分離奔月及環月軌道交會載人登月任務甘特圖
Fig.1 Gantt chart of manned lunar landing mission with crew and cargo trans-lunar separately and lunar orbit rendezvous

1.2 動力學模型

文獻[20]研究表明,地月轉移軌道偏差隨飛行時間增加而增益,著陸器地月轉移軌道飛行約5 d,應采用如下地心J2000坐標系中高精度軌道動力學模型:

(1)

式中:r為地心位置矢量;μE為地球引力常數;AN為N體引力攝動加速度,地月空間只需考慮日月攝動即可,相對位置通過JPL實驗室公布的DE405/LE405星歷求解;ANSE為地球非球形攝動,取WGS84引力場模型6×6階計算;ANSM為月球非球形攝動,取LP165P引力場模型6×6階計算;AR為太陽光壓,具體根據不同階段飛行器面值比和反射因子計算;AD為大氣阻力攝動加速度;AP為推進加速度。此處忽略木星等大行星的攝動、地球潮汐的攝動、地球扁率的間接攝動及相對論效應等微小量。

月心J2000坐標系中高精度軌道動力學模型與式(1)類似,即

(2)

式中:ρ為月心位置矢量;μM為月球引力常數;AN只需考慮日地攝動即可,相對位置通過JPL實驗室公布的DE405/LE405星歷求解;月心軌道段無大氣攝動項;其他項與式(1)相同。

2 窗口分層搜索策略

圖2 三層逐級搜索策略示意圖
Fig.2 Schematic of three-level cascade searching strategy

著陸器奔月窗口不受航天員生保系統影響,每個朔望月一般存在升段到達和降段到達兩次月窗口,直接采用循環結構搜索窗口,每次循環內部采用高精度軌道動力學模型求解奔月軌道存在兩個方面難題:一是循環步長大小產生漏解或計算量大矛盾;二是如果采用局部收斂算法求解高精度軌道動力學模型下結果,存在初值難猜測困難,如果采用不需初值的進化算法,必然千百倍增加計算量。本文提出的窗口三層逐級搜索策略如圖 2 所示。第一層采用地心二體解析公式,結合LEO傾角約束,在一個朔望月范圍內查找月窗口區間;第二層采用改進的雙二體假設解析模型,循環逆向快速求解窗口解集規律,獲取窗口初值;第三層采用SQP_Snopt算法和高精度軌道動力學模型,優化求解零窗口及奔月軌道精確值。

2.1 地心二體解析估計

(3)

(4)

(5)

(6)

2.2 改進雙二體假設解析模型打靶

相比于行星探測問題,地月距離較近、月-地質量比較大,雙二體模型求解地月轉移軌道與高精度模型相差數千公里,常用于月窗口搜索[11,19]和軌道參數特性分析等[9-10]。

(7)

tsph=tprl-ΔtMJ2

(8)

(9)

tTLI=tsph-ΔtEJ2

(10)

雙二體假設逆向求解奔月軌道如圖3所示。

圖3 雙二體假設逆向軌道求解示意圖
Fig.3 Schematic of double two-body reverse orbit joint

地月轉移總時長為

Δttotal=ΔtMJ2+ΔtEJ2

(11)

2.3 高精度軌道動力學模型優化求解

篩選改進雙二體假設解析模型解析打靶參數作為高精度軌道動力學模型下優化問題初值,采用SQP_Snopt算法迭代,優化變量為

(12)

優化變量上下界約束為

(13)

minJ=Q1J1+Q2J2+Q3J3

(14)

式中:Q1、Q2和Q3用來調整參數單位不同帶來的收斂不協調問題;J1和J2分別為TLI時刻近地距和軌道傾角優化目標;J3為控制飛行時長的真近點角優化目標。

(15)

式中:上標“*”表示瞄準值。

3 仿真分析

3.1 參數配置

以8 Apr 2025 18:35:00.000 UTCG時刻動力下降月面虹灣為例,設置月面活動3 d,參考文獻[16]生成該時刻100 km高環月圓軌道修正軌道六根數,如表1所示。表中:PeriRad為近地距;Ecc為偏心率;I為軌道傾角;RAAN為升交點赤經;ArgPeri為近拱點角距。

表1 環月圓軌道修正軌道六根數Table 1 Modified orbital elements of LLO

3.2 地心二體解析估計結果

采用式(2)將表 1 中LLO軌道參數高精度逆向積分至27 Mar 2025 18:35:00.000 UTCG時刻,再逆向積分至27 Feb 2025 18:35:00.000 UTCG時刻,得到LLO在一個朔望月內軌道面變化如圖4所示。可見高精度軌道動力學模型中,LLO軌道攝動演化規律不明顯,且變化范圍較大,采用高精度數值積分求解是必要的。

采用地心二體模型估算奔月軌道地心最小傾角如圖 5 所示。圖中:Moondeclination為月球赤緯。2025年白道面與赤道面相對夾角較大,約為28°36′,只有LEO傾角大于該值,才具備隨時奔月窗口。

以28.5° LEO傾角(參考西昌發射中心102°E/28.5°N)計算奔月軌道升交點赤經(RAAN)和近拱點角距(ArgPeri),如圖6所示,半長軸(SMajAX)和偏心率(Ecc)如圖 7 所示。

圖4 LLO軌道面攝動變化
Fig.4 Variation of LLO plane perturbation

圖5 地心二體模型求解的奔月窗口
Fig.5 TLI windows of geocentric two-body model

圖6 地心二體模型求解的升交點赤經和近拱點角距
Fig.6 TLI orbital RAAN and ArgPeri of geocentric two-body model

圖7 地心二體模型求解的半長軸和偏心率
Fig.7 TLI orbital SMajAX and E of geocentric two-body model

3.3 改進雙二體假設解析模型打靶結果

從圖8可以看出,一個朔望月存在兩個月窗口區間,窗口區間大小與LEO傾角正相關。地月轉移時長對朔望月周期不敏感。

3.4 高精度軌道動力學模型優化求解結果

考慮發射場指揮協調時間,圖 8(a) 中27 Mar 2025 18:35:00.000 UTCG至20 Mar 2025 17:50:00.000 UTCG月窗口不可用,設置搜索區間為11 Mar 2025 04:35:0.000~16 Mar 2025 04:35:00.000 UTCG,以2.3節所述優化模型和SQP_Snopt算法求解。得到精確奔月軌道時刻為9 Mar 2025 01:18:00.000 UTCG,時刻為3 Mar 2025 19:57:32.000 UTCG的一條軌道。對應月心J2000坐標系中近月點修正軌道六根數和地心J2000坐標系中TLI時刻修正軌道六根數,分別如表2和表3所示。

圖8 雙二體假設解析模型打靶TLI時刻地心近地距、地心軌道傾角和地月轉移時長
Fig.8 Geocentric PeriRad,inclination,and ΔT of TLI with double two-body model

參數與商業軟件STK校驗結果一致,空間軌跡如圖 10 所示。

仿真PC為Inter?CoreTMi5-2400 CPU @3.10 GHz 2.99 GB配置,采用VC++6.0代碼仿真,地心二體解析解估計和改進雙二體假設解析模型解析打靶仿真用時約40 s,高精度動力學逆向積分采用RKF7(8)積分器,優化計算時長與設置的迭代次數正相關,經過雙二體假設解析模型初值求解,算例設置迭代10次,Q1=1×10-6,Q2=Q3=1,計算用時190 s,目標函數收斂變化如圖 11 所示。

圖9 雙二體假設解析模型TLI時刻地心近地距、地心軌道傾角等高線圖
Fig.9 Geocentric PeriRad and geocentric inclination contour of TLI with double two-body model

表2 奔月軌道近月點修正軌道參數Table 2 Trans-lunar orbital modified perilune parameters

表3 奔月軌道近地點修正軌道參數Table 3 Trans-lunar orbital modified perigee parameters

圖10 著陸器奔月軌跡
Fig.10 LM trans-lunar trajectory

圖11 SQP_Snopt算法高精度動力學模型迭代目標優化
Fig.11 Optimization object iteration of SQP_Snopt algorithm with high precision dynamics model

4 結 論

本文構建了載人登月著陸器奔月窗口地心二體估計、改進雙二體假設解析模型和高精度軌道動力學優化求解模型,提出的三層逐級搜索方法能高效精確求解奔月零窗口。仿真結果表明:

1) 改進的雙二體簡化模型是有效的,提出的三層逐級搜索方法可以快速精確求解本應大量計算的著陸器奔月零窗口。

2) 將雙二體假設解析模型求解結果作為高精度軌道動力學模型進一步優化求解初值是必要的,SQP_Snopt算法在該策略中表現性能優異。

致 謝

感謝載人航天總體論證研究中心彭祺擘博士共同討論。

[1] COLAPRETE A, SCHULTZ P, HELDMANN J, et al. Detection of water in the LCROSS ejecta plume[J]. Science, 2010, 330(6003): 463-468.

[2] MITROFANOV I G, SANIN A B, BOYNTON W V, et al. Hydrogen mapping of the lunar south pole using the LRO neutron detector experiment LEND[J]. Science, 2010, 330(6003): 483-486.

[3] WU W R, LIU W W, QIAO D, et al. Investigation on the development of deep space exploration[J]. Science China: Technological Sciences, 2012, 55(4): 1086-1091.

[4] 李楨, 周建平, 程文科, 等. 環月軌道交會的奔月方案[J]. 國防科技大學學報, 2009, 31(1): 16-20.

LI Z, ZHOU J P, CHENG W K, et al. Investigation on lunar mission based on lunar orbit rendezvous [J]. Journal of National University of Defense Technology, 2009, 31(1): 16-20 (in Chinese).

[5] 周建平. 載人航天交會對接技術[J]. 載人航天, 2011, 17(2): 1-8.

ZHOU J P. Rendezvous and docking technology of manned space flight[J]. Manned Spaceflight, 2011, 17(2): 1-8 (in Chinese).

[6] 鄭愛武, 周建平. 載人登月軌道設計方法及其約束條件概述[J]. 載人航天, 2012, 18(1): 48-54.

ZHENG A W, ZHOU J P. A survey on trajectory design and constrains of manned lunar landing missions[J]. Manned Spaceflight, 2012, 18(1): 48-54 (in Chinese).

[7] 鄭愛武, 周建平. 直接再入大氣的月地返回窗口搜索策略[J]. 航空學報, 2014, 35(8): 2243-2250.

ZHENG A W, ZHOU J P. A search strategy of back windows for moon-to-earth trajectories directly re-turning to the earth[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(8): 2243-2250 (in Chinese).

[8] 周文艷, 楊維廉. 嫦娥二號衛星軌道設計[J]. 航天器工程, 2010, 19(5): 24-28.

ZHOU W Y, YANG W L. Orbit design for Chang’e-2 lunar orbiter[J]. Spacecraft Engineering, 2010, 19(5): 24-28 (in Chinese).

[9] PENG Q B, SHEN H X, LI H Y. Free return orbit de-sign and characteristics analysis for manned lunar mission[J]. Science China: Technological Sciences, 2011, 54(12): 3243-3250.

[10] SHEN H X, ZHOU J P, PENG Q B, et al. Point return orbit design and characteristics analysis for manned lunar mission[J]. Science China: Technological Sciences, 2012, 55(9): 2561-2569.

[11] LI J Y, GONG S P, WANG X. Launch window for manned moon-to-earth trajectories[J]. Aircraft Engineering and Aerospace Technology, 2012, 84(5): 344-356.

[12] 賀波勇, 沈紅新, 李海陽. 地月轉移軌道設計的改進微分校正方法[J]. 國防科技大學學報, 2014, 36(6): 60-64.

HE B Y, SHEN H X, LI H Y. An improved differential correction method for trans-lunar orbit design[J]. Journal of National University of Defense Technology, 2014, 36(6): 60-64 (in Chinese).

[13] 張洪禮, 羅欽欽, 韓潮, 等. UKF參數估計在三體Lambert問題中的應用[J]. 北京航空航天大學學報, 2015, 41(2): 228-233.

ZHANG H L, LUO Q Q, HAN C, et al. Application of UKF parameter estimation in the three-body Lambert problem[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(2): 228-233 (in Chinese).

[14] 羅欽欽, 韓潮. 包含引力輔助變軌的三體Lambert 問題求解算法[J]. 北京航空航天大學學報, 2013, 39(5): 679-682.

LUO Q Q, HAN C. Solution algorithm of the three body Lambert problem with gravity assist maneuver[J]. Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(5): 679-682 (in Chinese).

[15] PHILIP E G, MURRAY W, SAUNDERS M A. SNOPT: An SQP algorithm for large-scale constrained optimization[J]. SIAM Journal on Optimization, 2002, 12(4): 979-1006.

[16] CONDON G. Lunar orbit insertion targeting and associated outbound mission design for lunar sortie missions: AIAA-2007-6680[R]. Reston: AIAA, 2007.

[17] 賀波勇, 彭祺擘, 沈紅新, 等. 載人登月軌道月面可達區域分析[J]. 載人航天, 2014, 20(4): 290-295.

HE B Y, PENG Q B, SHEN H X, et al. Reachable region analysis of orbits for manned lunar landing mission[J]. Manned Spaceflight, 2014, 20(4): 290-295 (in Chinese).

[18] 沈紅新, 李海陽, 彭祺擘, 等. 探月飛行器定點返回軌跡特性分析[J]. 國防科技大學學報, 2011, 33(4): 6-11.

SHEN H X, LI H Y, PENG Q B, et al. Point return trajectory characteristics analysis for a lunar spacecraft[J]. Journal of National University of Defense Technology, 2011, 33(4): 6-11 (in Chinese).

[19] 賀波勇, 李海陽, 周建平. 載人登月繞月自由返回軌道與窗口精確快速設計[J]. 宇航學報, 2016, 37(5): 512-518.

HE B Y, LI H Y, ZHOU J P. Rapid design of circumlunar free-return high accuracy trajectory and trans-lunar window for manned lunar landing mission[J]. Journal of Astronautics, 2016, 37(5): 512-518 (in Chinese).

[20] 賀波勇. 載人登月轉移軌道偏差傳播分析與中途修正方法研究[D]. 長沙:國防科學技術大學,2013.

HE B Y. Analysis of transfer orbit deviation propagation and mid-course correction for manned lunar landing[D]. Changsha:National University of Defense Technology, 2013 (in Chinese).

[21] 郗曉寧, 曾國強, 任萱, 等. 月球探測器軌道設計[M]. 北京: 國防工業出版社, 2001: 53-58.

XI X N, ZENG G Q, REN X, et al. Orbit design of lunar probe[M]. Beijing: National Industry Press, 2001: 53-58 (in Chinese).

Lunarmoduletrans-lunarwindowsearchingapproachformannedlunarmission

HEBoyong,LIHaiyang*

CollegeofAerospaceScienceandEngineering,NationalUniversityofDefenseTechnology,Changsha410073,China

Anaccurateandrapidapproachtolunarmodule(LM)trans-lunarinjection(TLI)windowandorbitparameterdesignisderivedformannedlunarlandingmissionbasedonlowlunarorbitcoplanarrendezvousanddocking(RVD).Thepatternthatastronautsareseparatedfromcargoisadopted.LMarrivesatthelunarcoplanarRVDorbitearlierthanthecrewexplorationvehicle(CEV).LMcompletesthebrakingontheLLORVDplanewithonlyoneimpulseatperilune.Astrategyofthree-levelcascadesearchingforLMTLIwindowsisproposed.ThemonthwindowandorbitalinitialparametersofTLIareobtainedinthefirstlevel.Furtherinitialsolutionformonthwindowandorbitalparametersareobtainedwithanimproveddoubletwo-bodyanalyticmodelinthesecondlevel.ThehighprecisiondynamicssolutionforzerowindowandorbitalparametersareoptimizedbySQP_Snoptalgorithm.Simulationresultsshowthattherapidandaccuratethree-levelcascadesearchingapproachforLMTLIwindowsandorbitalparametersiseffectiveformannedlunarlandingmission.TheinfluentialfactorsandfundamentalpatternsofTLIwindowarerevealedtosupply

toChina’smannedlunarlandingmissioninthefuture.

mannedlunarlanding;lunarmodule(LM);lowlunarorbitrendezvous;trans-lunarwindow;cascadesearching

2016-06-30;Revised2016-08-17;Accepted2016-09-26;Publishedonline2016-10-091007

URL:www.cnki.net/kcms/detail/11.1929.V.20161009.1007.008.html

s:NationalNaturalScienceFoundationofChina(11372345,11402295);NationalBasicResearchProgramofChina(2013CB733100)

2016-06-30;退修日期2016-08-17;錄用日期2016-09-26; < class="emphasis_bold">網絡出版時間

時間:2016-10-091007

www.cnki.net/kcms/detail/11.1929.V.20161009.1007.008.html

國家自然科學基金 (11372345,11402295); 國家“973”計劃 (2013CB733100)

.E-maillihaiyang@nudt.edu.cn

賀波勇, 李海陽. 載人登月著陸器奔月窗口搜索方法J. 航空學報,2017,38(4):320574.HEBY,LIHY.Lunarmoduletrans-lunarwindowsearchingapproachformannedlunarmissionJ.ActaAeronauticaetAstronauticaSinica,2017,38(4):320574.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0261

V412.4+1

A

1000-6893(2017)04-320574-09

(責任編輯: 張玉)

*Correspondingauthor.E-maillihaiyang@nudt.edu.cn

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产精品一区二区国产主播| 午夜精品影院| 狠狠躁天天躁夜夜躁婷婷| 国产理论精品| 亚洲精品午夜天堂网页| 亚洲成a人在线播放www| 99国产精品国产高清一区二区| 亚洲第一黄色网| 55夜色66夜色国产精品视频| 欧美亚洲国产一区| 91九色最新地址| 久久综合丝袜长腿丝袜| 婷婷色在线视频| 91视频青青草| 国产精品午夜福利麻豆| 欧美自慰一级看片免费| 久久久久国产精品嫩草影院| 无码AV动漫| 野花国产精品入口| 欧美另类一区| 国产剧情一区二区| 五月丁香在线视频| 国外欧美一区另类中文字幕| 国产对白刺激真实精品91| 国产精品一区在线麻豆| 欧美成人手机在线观看网址| 欧美日韩高清在线| 强奷白丝美女在线观看| 欧美亚洲国产日韩电影在线| 久久综合一个色综合网| 欧美一级夜夜爽www| 亚洲美女操| 欧美一区福利| 秋霞一区二区三区| 亚洲无线国产观看| 午夜三级在线| 亚洲成aⅴ人在线观看| 欧美日韩中文字幕在线| 夜夜操天天摸| 91精品啪在线观看国产91| 国产成人无码AV在线播放动漫| 五月天久久综合| 婷婷99视频精品全部在线观看| 中文字幕精品一区二区三区视频| 福利小视频在线播放| 欧美中文字幕一区二区三区| 国产亚洲精品yxsp| 国产精品一老牛影视频| 国产成人福利在线视老湿机| 2019年国产精品自拍不卡| 人与鲁专区| 久久天天躁狠狠躁夜夜2020一| 亚洲av无码专区久久蜜芽| 国产亚洲精品97AA片在线播放| 国产精品亚洲综合久久小说| 91网站国产| 女人一级毛片| 日本91视频| 久久亚洲精少妇毛片午夜无码| 亚洲日产2021三区在线| 天天综合网色| 久久熟女AV| 精品一区二区三区自慰喷水| 亚洲av无码牛牛影视在线二区| 亚洲天堂网在线视频| 黄色网站不卡无码| 国产一区二区三区视频| 激情综合网激情综合| 极品私人尤物在线精品首页| 一边摸一边做爽的视频17国产| 亚洲手机在线| 99re在线观看视频| 91色在线观看| 欧美精品在线看| 欧美色丁香| www.精品视频| 欧美人在线一区二区三区| 精品视频第一页| 国产白丝av| 日本高清视频在线www色| 亚洲va视频| 国产成年无码AⅤ片在线|