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

一種多級全固體運載火箭上升段自主制導方法

2019-02-15 08:51:42李新國
宇航學報 2019年1期
關鍵詞:發動機

張 遷,許 志,李新國

(1. 西北工業大學航天學院,西安 710072;2. 陜西省空天飛行器設計技術重點實驗室,西安 710072)

0 引 言

近年來小衛星由傳統的搭載開始向商業購買發射服務轉變,呈現個性化、定制化、眾籌化的特點,而小型固體運載火箭憑借快速機動發射、可靠性高、成本低等特點業已成為發射任務的首要選擇。2013年“Epsilon”火箭的成功發射,在快速響應航天運載器領域產生了很大震動,其標準型采用了三級全固體火箭方案[1]。由于免去了液體火箭在臨射前加注燃料等流程,最短發射準備時間縮短至24 h以內,在自然災害、突發事件的應急發射任務中具有顯著優勢,同時在商業航天領域同樣有著強大的競爭力[2]。然而,相對于液體火箭發動機,固體火箭為了提高質量比及可靠性取消了推力終止機構,導致只能采取燃料耗盡關機而不能進行制導關機;此外,固體發動機具有工作時間短、推力大等特點,必須采用“助推-滑行-助推”飛行模式才能保證對不同發射任務的適應性。上述特點給全固體火箭的制導技術帶來了新的困難與挑戰。

國內外學者針對上升段制導算法的研究主要集中在液體運載火箭的最優閉環制導問題。Brown等[3]針對大氣層外飛行器的在線制導給出了線性模型最優解析解。在此基礎上,文獻[4]設計了一種迭代制導算法在線迭代求解最優解析解的待定參數;文獻[5]對閉環最優制導進行梳理,總結了多種任務約束下迭代制導方法的最新研究成果。文獻[6-8]等根據變分法最優控制原理推導了燃料最省的上升段軌跡最優解,并分別采用了多重打靶法、有限差分法、配點法等數值算法進行求解。文獻[9-10]開展了固體火箭上升段彈道快速設計方法研究,設計了基于Gauss偽譜法的軌跡優化策略使軌跡優化效率進一步提高。對于液體火箭或具有液體上面級的固體火箭,上述算法能夠在滿足終端約束的同時保證飛行軌跡的最優性;然而,針對耗盡關機的全固體運載火箭,必須考慮剩余能量管控與發動機大散差抑制等問題,這與以燃料最省的液體運載火箭制導方法有顯著的差別。因此,在快速響應需求下全固體運載火箭制導必須解決具有固定總沖約束的兩點邊值問題。

而固體火箭制導算法的研究主要是對經典的閉路制導[11-13]進行改進。針對閉路制導魯棒性問題,文獻[11]考慮發動機性能參數散布特性,提出需用速度增益曲面的概念并設計了基于該概念的大氣層外導引方法。針對以軌道根數為終端約束的制導問題,文獻[12]根據實時飛行狀態與軌道根數之間的關系推導出需要速度求解方程,并應用于發射近圓軌道的閉路制導控制中。文獻[13]考慮了地球非球形J2攝動項對彈道的影響,改進了需要速度求解算法以修正地球引力攝動在長時段飛行過程中引起的求解偏差。由于未考慮耗盡關機方式下火箭的速度管控問題,上述閉路制導算法無法直接應用于具有軌道能量約束的太陽同步軌道(SSO)發射任務。針對耗盡關機能量管理問題,文獻[14-16]等通過設計不同的能量管理曲線來實現對耗盡關機速度矢量的管控。但是,上述算法由于以需要速度矢量為導引量而忽略了制導過程中位置矢量變化對終端約束產生的影響,故針對終端多約束條件的入軌任務其制導精度較差,魯棒性無法保證。

為滿足終端多軌道根數約束,文獻[17-20]根據軌道動量矩守恒定律推導出具有速度矢量及位置矢量約束的定點制導算法。本文在PA[20]制導算法的基礎上設計了一種適應于多子級“助推-滑行-助推”模式的制導算法,通過具有軌道能量約束的滑行過渡軌道將多級推進段制導統一起來,根據當前飛行狀態與終端軌道根數約束推導出滿足發動機耗盡關機條件的過渡軌道求解方程,進而確定出軌道根數約束與火箭點火時間及推力矢量方向的數學關系。

1 固體火箭上升段制導問題

1.1 質心運動模型

在發射慣性坐標系下運載火箭質心運動方程[20]表示為:

(1)

T為發動機推力,軸向和法向氣動力分別為A和N,其表達式為:

A=qSCA,N=qSCN

xb和yb分別表示彈體軸向和法向。固體運載火箭發動機的推力大小與秒流量不可調節,式(1)中實際上需要確定出彈體方向xb的控制指令,使運載火箭以耗盡關機的方式滿足入軌約束條件。當運載火箭進入真空環境,大氣壓強P(h)產生的推力及火箭所受的氣動力予以忽略,運載火箭所受的氣動載荷及氣動角等過程約束也隨之消失,為自主制導算法的設計提供了更大的空間。由于大氣層內飛行軌跡的偏差可以通過大氣層外自主制導算法進行修正,本文重點研究大氣層外多約束自主制導算法。

1.2 實際推力非線性問題

固體發動機能夠提供的視速度模量WM隨推進劑質量ms變化的表達式為:

(2)

在耗盡關機模式下,當推進劑質量ms確定則視速度模量WM為常值。視位置模量RM隨推進劑質量ms變化的表達式為:

(3)

發動機持續推進過程的形心時間Tc表達式為:

(4)

然而,發動機在燃燒過程中由于外界壓強、溫度等環境因素以及燃燒過程化學反應程度的影響,實際性能參數與工況條件下相比偏差較大,非線性程度較液體發動機更為明顯,如圖1所示。

圖1 固體發動機性能參數曲線在高低溫條件下散布示意圖Fig.1 Profile of the performance distribution with solid motor thrust-time diagrammatic drawing

當發動機裝藥量ms一定時,圖1中不同條件下曲線所圍的面積是相同的,式(2)依然成立,即視速度模量為常值,但視位置模量隨著秒流量的散布而產生偏差。由于火箭在實際飛行過程中發動機真實燃燒曲線無法提前獲知,發動機模型的不確定性及非線性使制導律的設計更為困難,最優化制導算法通常未考慮實際非線性推力散布模型或者認為發動機的模型提前已知,因此無法直接應用于實際工程問題。針對上述發動機特性,要求制導律在滿足約束條件下盡可能減小姿態角變化率以增加終端制導精度。

1.3 耗盡關機多約束問題

運載火箭的終端約束為衛星載荷的軌道要素,通常軌道根數約束為半長軸a、偏心率矢量e及軌道傾角i等,表達式如下:

(5)

固體運載火箭入軌任務約束為軌道要素方程式(5),由于終端軌道根數約束條件的不同,運載火箭受到約束的狀態量也隨之改變。此外,由于固體火箭發動機以燃料耗盡的方式關機,推進段發動機產生的總速度增量和總位置增量為定值,滿足式(2)及式(3)約束。

固體運載火箭由于在大氣層外通常采取“助推-滑行-助推”飛行模式,即在滑行階段確定發動機點火時間,在推進階段確定制導指令,使式(1)滿足終端約束。因此,針對全固體運載火箭,制導算法本質上是求解具有固定總沖約束的兩點邊值問題,而其算法的核心則是:

(1) 推導出推力矢量、點火時間與終端狀態約束之間的理論關系;

(2) 確定能量匹配的滑行過渡軌道以保證推進段耗盡關機約束要求。

2 PA制導算法

2.1 PA理論基礎

固體火箭發動機推力脈沖大,工作時間短,可以快速產生瞬時速度沖量,此時火箭彈體方向沿著速度沖量方向[11-13]。實際上,發動機需要經過額定工作時間才能完全耗盡產生相應的速度沖量,PA制導算法主要研究“瞬時脈沖矢量”與“定向持續推力矢量”在實現入軌任務時的理論關系[17-20]。因此,假設發動機持續推進過程中推力矢量方向始終沿著速度沖量方向Γ。詳細理論推導過程參考文獻[20]。

當運載火箭在真空環境中以定向的推力矢量方向飛行,其動力學方程為:

(6)

式中,由引力產生的狀態分量表達式為:

式(6)所描述的運載火箭運動方程可分解為:沿原來的開普勒軌道繼續滑行及沿定向推力的作用產生位移矢量和速度矢量,表達式如下:

(7)

那么,持續推力過程引起運載火箭動量矩的變化為:

ΔH=mfrorb, f×vorb, f-m0r0×v0

(8)

將式(7)代入式(8)得到:

ΔH=mf(rsub, f×vsub, f+rsub, f×WMΓ+

RMΓ×vsub, f+RMΓ×WMΓ)-m0r0×v0

(9)

化簡式(9)并整理得:

ΔH=mfrsub, f×vsub, f-m0r0×v0+

mf[rsub, f-(RM/WM)vsub, f]×(WMΓ)

(10)

在過渡軌道面和目標軌道面確定出軌道面交線Rimp,軌道相關參數的矢量關系如圖2所示。

圖2 PA制導算法運動分解圖Fig.2 Motion analysis of pointing guidance method

令等效位置矢量rimp為:

rimp=rsub, f-(RM/WM)vsub, f

(11)

將式(11)代入式(10)得:

ΔH=mfrimp×(vsub,imp+WMΓ)-m0rimp×vimp

(12)

式(12)表明:“定向持續推進過程”對軌跡的改變與在等效脈沖點Pimp處施加“瞬時脈沖矢量”對軌跡的影響等效。

2.2 PA制導指令解算

運載火箭的狀態矢量(r,v) 通常建立在地面發射坐標系,通過地心慣性赤道坐標系與發射系之間的變換矩陣得到火箭的慣性狀態矢量,從而得到相應的軌道根數。在軌道面交點Pimp處,速度矢量關系如圖3所示,通過滑行軌道參數求得:

圖3 軌道交點處速度矢量關系Fig.3 Vector diagram of the orbit intersection point

(13)

同樣地,根據目標軌道參數得到:

(14)

根據式(13)和式(14),求解出軌道面交點Pimp處的真近點角、切向速度及法向速度。同時根據軌道傾角i和Pimp的經度、緯度,求解出速度矢量的方位角Asub和Aorb。因此,在兩軌道交點Pimp處速度矢量表示如下:

(15)

進而得到PA制導所需的速度矢量大小和方向:

(16)

軌道面交點Pimp根據終端軌道根數約束而求得,則火箭從當前點飛行至交點Pimp處的時間為:

(17)

timp-sub, f=timp+RM/WM

(18)

根據式(4)、式(17)和式(18)得運載火箭點火時間為:

ta-ig=ta-imp-Tc

(19)

綜上所述,以定向推力矢量的方式進行制導時,終端軌道位置矢量約束通過滑行軌道的點火時間式(19)來滿足,速度矢量約束通過發動機定向推力矢量式(16)來滿足,從而滿足終端多軌道根數約束。

3 多子級改進PA算法

3.1 PA理論拓展

對于耗盡關機的固體運載火箭,由于無法制導關機,火箭對速度、位置的管控能力受限。為滿足終端多軌道要素約束條件,單子級火箭采取附加姿態機動的方式耗散多余能量,導致姿態角速率變化劇烈;多子級火箭可以通過“助推-滑行-助推”的模式,使火箭的動能與勢能轉換,通過最佳的滑行軌道來匹配各級火箭發動機的總沖量,能量管控的實現機理如圖4所示。

圖4 多子級改進PA算法能(MPA)量匹配過程示意圖Fig.4 Profiles of the energy matching process by the MPA

在式(12)基礎上,多子級火箭持續推力過程引起的軌道動量矩變化為:

ΔH=mnrn×vn-…-m1r1×v1-m0r0×v0

(20)

根據PA理論,將式(12)代入式(20)得:

ΔH=mnrn,imp×(vn,imp+WM,nΓn)-…-

m1r1,imp×(v1,imp+WM,1Γ1)-m0r1,imp×v1,imp

(21)

各級火箭發動機推進過程中動量矩的變化通過脈沖點ri,imp(i∈[1,n]) 聯系起來。為滿足耗盡關機條件,由式(21)、式(16)及式(2)得:

Wi,PA=WM,i,i=1,2,…,n

(22)

根據終端軌道根數約束式(5),通過改變定向推力矢量方向Γi(i∈[1,n])來達到耗盡關機終端多約束條件。火箭各階段的點火時間根據式(17)得:

(23)

3.2 改進算法降階求解

以PA理論為基礎,滑行過渡軌道將多子級耗盡關機問題統一起來,根據滑行軌道參數解算推力矢量方向Γ。以四子級全固體運載火箭發射SSO軌道為例,求解最佳能量匹配滑行軌道參數。

根據終端軌道根數條件及SSO軌道偏心率eorb為零得:

(24)

根據開普勒軌道的性質,火箭處于滑行軌道階段,軌道動量矩和軌道能量守恒,因此當前軌道參數與運載火箭狀態量滿足:

(25)

同樣地,滑行過渡軌道參數與運載火箭狀態量滿足:

(26)

其中,v0⊥,v1,orb⊥,v2,sub⊥分別表示在點P0,P1,imp,P2,imp處的徑向速度分量。在軌道交點Pi,imp處速度矢量方向如圖3所示,相應的表達式為:

(27)

由于軌道傾角的調整發生在P1,imp處完成,火箭的后續飛行將處于軌道面內,將終端約束式(24)代入式(26),并根據矢量方向關系式(27)得:

(28)

同樣地,將滑行過渡軌道約束式(28)代入式(25),并根據矢量關系式(27)得:

(29)

根據運載火箭導航解算的狀態向量(r,v) 及等式方程組式(25)、式(26)、式(28)和式(29),解算滿足耗盡關機多約束條件下的推力矢量方向Γ。上述方程組為超越方程無法直接獲得解析解,為保證算法的有效性,以v1,orb為迭代變量,將非線性方程組的求解問題轉化為一維求根問題,具體流程如下:

(1) 確定迭代變量v1,orb區間[v1,sub,v1,sub+WM,1]:由于火箭處于上升段有v1,orb≥v1,sub,根據速度矢量得v1,orb≥v1,sub+WM,1。

(2) 解算推力矢量方向Γ2,從而得到切向速度v2,sub⊥;根據軌道能量方程式(26)和約束式(28),得:

(30)

(31)

因此,更新后的速度為:

cos(?1,sub-Γ1,φ)

(32)

(5) 根據v1,orb和v2,sub⊥解算最佳能量匹配滑行軌道參數:

(33)

運載火箭在實際飛行過程中,由于受到干擾及不確定性的影響或發生非災難性故障導致原定飛行任務無法現實。在這種特殊情況下,運載火箭由于后續運載能力的不足,使原方程的收斂解處于流程(1)的迭代區間之外,導致迭代算法失效。針對上述問題,制導算法需要根據當前運載火箭的飛行條件在線調整飛行任務約束以保證火箭的飛行安全,飛行任務的終端軌道根數調整過程如圖5所示。

圖5 飛行任務調整示意圖Fig.5 The adjustment process of trajectory geometry

根據軌道根數表達式(5)可知,軌道半長軸是軌道能量的函數,因此,通過降低軌道半長軸使原方程的收斂解處于流程(a)的迭代區間內,以保證運載火箭盡可能地實現原定飛行任務。對于一維求根問題,當方程的解在迭代區間內時,黃金分割法在迭代10~20步以內可獲得高精度解,根據式(33)計算滿足能量約束的滑行過渡軌道參數,然后根據式(16)解算導航坐標系內的制導指令,并根據式(19)計算火箭點火時間。隨著火箭持續飛行,多級問題將轉換為單級問題,根據式(16)和式(19)在線解算制導指令,直至火箭進入目標軌道。

4 仿真校驗與性能評估

飛行任務以太陽同步軌道(SSO)為例,考慮地球自旋角速度及引力J2項攝動。本文研究對象采用四子級全固體運載火箭,其制導算法流程如圖6所示,制導方案如下所述。

1) 大氣層內采用開環方案飛行,二級耗盡關機后,運載火箭位于大氣層外飛行軌跡的偏差以導航輸入的方式提供給自主制導算法。

2) 運載火箭第三子級采用多級改進PA算法,在線解算出最佳能量匹配軌道,通過調整滑行過渡軌道的能量與動量矩,使發動機總能量與飛行任務相匹配,以彌補發動機采用耗盡關機時運載火箭對速度矢量管控的不足。

3) 由于第三子級已經考慮了發動機多余能量的消耗問題,運載火箭第四子級只需通過PA制導算法確定實際點火時間及推力矢量方向,即可使運載火箭以耗盡關機的方式進入目標軌道。

4) 第四子級耗盡關機時,火箭基本滿足入軌參數約束,末修級由于燃料的限制,主要小幅度修正速度偏差。

圖6 制導方案流程圖Fig.6 Flow chart of the guidance scheme

制導算法所需要的諸元參數按標稱值預先裝訂,在仿真模型中偏差散布及不確定性對飛行軌跡的影響以導航輸入的方式供制導算法在線解算制導指令,各項隨機偏差由每次打靶隨機產生,且服從正態分布,其中發動機模型采用實際內彈道推力曲線,并通過高低溫的形式表征發動機性能散差,模型不確定性及散差分布見表1。

表1 蒙特卡洛仿真散差配置表Table 1 Dispersions in Monte Carlo simulations

4.1 多軌道任務自適應性分析

為驗證制導算法對發射任務的適應性,考慮700 km SSO軌道最大載荷、空載荷兩種典型情況,并以最大載荷方案設置終端軌道高度為300 km、700 km兩種常用衛星軌道,終端狀態約束見表2。根據所設置的三種典型發射任務條件,各方案分別進行1000次蒙特卡洛仿真,其散差配置見表1。其中,火箭軌道地心距變量以終端地心距約束值進行單位化,絕對速度變量以終端入軌速度約束值進行單位化,三種發射任務條件下的蒙特卡洛仿真曲線簇繪制結果如圖7所示。

表2 運載火箭發射任務配置表Table 2 Missions and terminal state constraints in Monte Carlo

圖7(b)絕對速度變化曲線簇中出現飛行過程絕對速度超過目標速度的情況,因為300 km滿載與700 km空載分別對應軌道高度的變化與有效載荷質量的改變,使運載火箭均存在發動機能量過剩的情況,通過在線自主制導算法使狀態量均達到目標值,表明制導算法對不同任務的適應性,且針對固體火箭發動機不同程度的能量過剩情況均能達到耗盡關機條件。

圖7(a)~圖7(d)中各狀態量曲線隨時間推移均達到目標值,表明了制導算法的終端約束能力;

圖7 蒙特卡洛仿真狀態矢量曲線簇Fig.7 Profiles of the states run for Monte Carlo simulation

雖然偏差干擾源使運載火箭飛行軌跡改變,但應用在線自主制導算法后狀態量通過不同的曲線形式均達到目標值,表明了制導算法對終端狀態的高精度約束能力,其中曲線簇的寬度代表了制導算法對各干擾源的魯棒性。

4.2 制導精度與魯棒性分析

運載火箭以軌道根數為終端約束條件,飛行狀態矢量(r,v) 與軌道根數滿足相互對應關系式,但不同的軌道根數對狀態矢量的敏感度不同。軌道半長軸受速度大小的影響比較明顯,軌道偏心率與軌道傾角對速度矢量方向更為敏感。第4.1節仿真結果校驗了PA制導算法對終端狀態的約束能力,而針對多子級間軌道能量匹配與耗盡關機問題,通過運載火箭飛行過程中軌道能量與動量矩的變化,來驗證多級PA改進算法的自主性、制導精度及魯棒性。軌道能量與動量矩的變化如圖8所示。

圖8 多級PA自主能量匹配仿真曲線簇Fig.8 Profiles of the energy matching process by the multistage pointing algorithm

針對不同的發射任務,通過多級PA改進算法在線解算滑行過渡軌道根數約束條件,使發動機產生的速度矢量增量與實現軌道任務所需的速度矢量相匹配,來滿足固體火箭耗盡關機條件。在不同初始條件下圖8(a)軌道能量變化和圖8(b)動量矩變化均達到終端約束值,展現出制導算法的自主性,同時也體現出算法對不同干擾的魯棒性。

由圖9可知,終端偏差比較集中且散布較小,但依然存在野點,半長軸偏差的最大值為1520.9 m,軌道傾角偏差的最大值為1.4×10-2(°)。主要原因為:一方面,大氣層內未加橫法向導引的開環飛行,受氣動及風干擾影響導致飛行軌跡偏差較大;另一方面,發動機性能參數的偏差既不能提前獲知也無法在線測量,導致耗盡關機點必然存在狀態偏差。

3000次蒙特卡洛仿真的終端軌道要素偏差的統計結果見表3,不同發射任務條件下仿真曲線簇的終端入軌根數偏差如圖9所示。根據表3可得,半長軸偏差的數學期望值小于百米級,標準差小于200 m;偏心率的數學期望值與標準差達到10-4量級;軌道傾角的數學期望值達到10-3(°)量級,標準差達到10-2(°)量級。末修級燃料的平均值為6.5757 kg,最大值為14.761 kg,標準差小于1.7213 kg。數據結果表明,本文所提制導算法對制導關機的液體末修發動機依賴不強,對耗盡關機的固體發動機依然具有高精度終端約束能力。

圖9 終端軌道要素偏差散布Fig.9 Drogue deployment accuracy

5 結 論

1) PA制導算法根據軌道動量矩矢量方程并引入了速度矢量及位置矢量約束條件,推導出點火時間令與推力矢量求解方程,具有終端多軌道要素約束能力。

表3 蒙特卡洛仿真統計結果Table 3 Monte Carlo simulations statistic results

2) 多級改進PA算法,通過具有軌道能量約束的滑行過渡軌道建立多級能量匹配方程組,針對不同的發射任務具有在線解算制導指令并滿足耗盡關機多終端約束的能力。

3) 在載荷最大條件下,運載火箭幾乎無剩余能量,因此消耗末修級燃料較多;而在剩余能量過多的條件下,干擾源的影響使耗盡關機點存在狀態偏差并通過滑行軌道進一步放大,使制導算法的魯棒性降低。

4) 整體制導方案具有很高的制導精度,對所配置的干擾源具有很強的魯棒性,對不同的發射任務要求具有很強的自主性。

猜你喜歡
發動機
元征X-431實測:奔馳發動機編程
2015款寶馬525Li行駛中發動機熄火
2012年奔馳S600發動機故障燈偶爾點亮
發動機空中起動包線擴展試飛組織與實施
奔馳E200車發動機故障燈常亮
奔馳E260冷車時發動機抖動
新一代MTU2000發動機系列
2013年車用發動機排放控制回顧(下)
VM Motori公司新型R750發動機系列
發動機的怠速停止技術i-stop
主站蜘蛛池模板: 亚洲一欧洲中文字幕在线| 天天综合色网| YW尤物AV无码国产在线观看| 色国产视频| 久久一本日韩精品中文字幕屁孩| 91小视频版在线观看www| 日本免费福利视频| 亚洲综合一区国产精品| 一级做a爰片久久毛片毛片| 国产激情无码一区二区免费| 强奷白丝美女在线观看| 久久久久亚洲Av片无码观看| 久久永久视频| 精品视频在线观看你懂的一区| 国产一级二级三级毛片| 99999久久久久久亚洲| 亚洲h视频在线| 免费一级成人毛片| 一本久道久久综合多人| 久久特级毛片| 欧美国产日产一区二区| 99久久国产自偷自偷免费一区| 国产午夜一级毛片| 久久久受www免费人成| 亚洲无码高清一区| 亚洲国产综合精品中文第一| 日韩精品一区二区深田咏美| 99久久精品美女高潮喷水| 亚洲专区一区二区在线观看| 亚洲天堂久久久| 熟妇人妻无乱码中文字幕真矢织江| 久草国产在线观看| 四虎综合网| 在线免费观看a视频| 欧美精品高清| 久久婷婷五月综合色一区二区| 亚洲AV成人一区国产精品| 国产一级精品毛片基地| 深夜福利视频一区二区| 在线播放91| 在线99视频| 国产97色在线| 亚洲毛片一级带毛片基地| 精品三级网站| 国产精品爆乳99久久| 欧美伊人色综合久久天天| 九九九久久国产精品| 人妻熟妇日韩AV在线播放| 久久香蕉国产线看观看精品蕉| 伊人久久大线影院首页| 日本在线亚洲| 免费一级毛片不卡在线播放 | 97狠狠操| 538精品在线观看| 91成人免费观看在线观看| 在线观看免费黄色网址| 亚洲综合香蕉| 日韩欧美91| 天天色综网| 国产黑丝一区| 久久青青草原亚洲av无码| 日本久久免费| 婷婷激情亚洲| 国产欧美网站| 亚洲欧美人成人让影院| 国产偷倩视频| 波多野结衣视频一区二区 | 青青操国产视频| 国产sm重味一区二区三区| 亚洲综合欧美在线一区在线播放| 91麻豆国产视频| 2021最新国产精品网站| 97在线免费| 狠狠做深爱婷婷综合一区| 在线亚洲小视频| 中文字幕精品一区二区三区视频| 国产免费久久精品99re丫丫一| 在线观看91香蕉国产免费| 91无码视频在线观看| 亚洲第一成年网| 欧美亚洲一区二区三区导航| 啪啪免费视频一区二区|