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

超低軌衛星氣動舵機輔助姿態控制方法設計

2023-03-31 07:42:58王濤焦洪臣劉杰陳樂宇張迎春
北京航空航天大學學報 2023年3期

王濤,焦洪臣,劉杰,陳樂宇,張迎春

(1.哈爾濱工業大學 航天學院,哈爾濱150001;2.中國空間技術研究院 遙感衛星總體部,北京100094;3.北京航空航天大學 宇航學院,北京 100191)

超低軌衛星存在傳統衛星不可比擬的優勢。超低軌衛星軌道高度低,能夠快速到達預定軌道,展開工作;搭載的空間相機或雷達等設備指標無需太高,就能夠達到傳統衛星的效果,甚至更好,有效載荷這一項就能夠大幅度減少衛星成本;由于體積、質量和成本相對傳統遙感衛星可大幅縮減,超低軌衛星的部署和發射方式極為靈活,并且可通過多星組網的方式實現對特定區域優異的時間、空間分辨率[1-4]。

在早期的研究中,氣動力和力矩通常被視為干擾項,需要進行控制或低效,在實際的工程項目中也是如此應對。美國研發的NanoEye 衛星是一種用于對地觀測的圖像衛星,其運行的軌道高度范圍為200 ~300 km,外形設計考慮到超低軌衛星所受氣動力的影響,主結構為圓柱形,而在前端裝有2 塊太陽能電池板,構成一個楔子的形狀,減小作用在衛星上的氣動力[5]。俄羅斯的Yantar 系列及Orlets系列衛星均運行在超低軌道上,近地點高度約為200 km,遠地點高度約為300 km。此外,俄羅斯的部分“琥珀”系列衛星也運行在超低軌道上,為了抵抗氣動力導致的軌道衰減,第四代中的“琥珀-2K”衛星定期提升軌道,從而使近地點高度保持在170~180 km,遠地點高度保持在320~350 km[6]。受氣動阻力影響,這些衛星具有較強的對地觀測能力,但工作壽命較短。

隨著航天技術的發展,研究人員開始進行利用氣動力輔助超低軌衛星姿態控制方面的研究[7-11]。Kumar 等[12-13]研究了氣動被動穩定控制,對氣動穩定和磁阻尼穩定進行仿真,證明了氣動被動穩定控制的可行性。Psiaki[14-15]研究微納衛星的控制時,利用氣動力矩和磁力矩進行被動穩定控制,仿真結果表明,氣動力矩和磁力矩控制系統對于微納衛星姿態控制效果良好。Guettler[16]在2007 年提出利用衛星受到的氣動力矩實現衛星姿態主動控制,通過線性模型對衛星姿態運動建模仿真,證明了在較低軌道上才具有通過氣動力矩進行姿態控制的可行性。此外,研究人員發現低軌下的氣動力還可用于衛星動量輪角動量卸載[17-19]。歐洲航天局于2009 年3 月17 日發射GOCE 衛星,該衛星發射高度為275 km,在軌運行高度為250 ~260 km,主要用于重力場和海洋環流探測。衛星機身呈八邊形棱柱體,長約5.3 m,橫截面為1.1 m2。在衛星的后部安裝有2 對尾翼,當衛星姿態有偏移時,尾翼可以產生氣動力矩以修正姿態,但由于超低軌道空氣密度較低,這種被動控制方式調整緩慢,只是被用于輔助主動控制[20]。日本于2017 年12 月23 日發射SLATS 衛星,該衛星首先保持271.5 km 的軌道高度,隨后利用氣動力制動,最終在167.4 km 高度實現短期運行[21]。

針對滿足高分辨率、高重訪、全天時和快速響應需求的超低軌遙感衛星,研究滿足其空間環境和任務需求差異化特點的姿態控制策略。不同于以往將空間大氣環境對衛星的影響視為攝動干擾,通過補償技術對其進行修正,本文引入航空空氣動力學相關理論對控制手段和方法進行必要的補充和改進,在衛星姿態控制策略方面進行創新性研究,以氣動舵機作為輔助手段,有效降低衛星傳統控制機構壓力,實現超低軌環境下衛星姿態的高精度保持。

1 超低軌衛星氣動力和氣動力矩特性

1.1 大氣密度模型

本文采用航天領域常用的大氣模型——COSPAR 國際標準大氣模型(CIRA)來確定大氣密度,這是國際上搜集來的,通過測量成千上萬顆衛星由于大氣阻力造成的軌道衰減的效應的數據,最終加以確定,高度在500 km 以下的大氣密度可以采用指數衰減模型進行計算。指數衰減模型為

后續仿真分析中,大氣密度模型參照上述模型在程序中預先設定,依據衛星軌道高度實時計算當前大氣密度,并應用到氣動計算模型中,以求得較為精確的氣動力和力矩。

1.2 氣動計算模型

超低軌衛星運行環境屬于高層大氣,氣體稀薄,分子平均自由程遠大于衛星特征長度,因此,對于衛星受到的氣動力與氣動力矩模型按照自由分子流的理論進行分析。此時,衛星受到的氣動力和力矩依賴于多種因素,如大氣密度、衛星構型、表面材料特性、來流速度及矢量方向等。

在稀薄氣體動力學的3 個基本假設的前提下,基于氣體分子與表面相互作用的麥克斯韋模型,可得到表面壓力與剪切應力如下[22-23]:(2kT∞/m)1/2

式中:ρ∞為來流密度;V∞為來流速度;Tw/T∞為壁面溫度和來流溫度比;S 為分子速度比,定義為S= V∞/,k 為玻爾茲曼常數,m 為單個分子質量;σn和στ分別為法向和切向動量適應系數;θc為來流和表面的夾角;erf(Ssin θ)為誤差函數,定義為

氣動力計算采用當地化方法,基本原理與稀薄氣體動力學相同,只是定義來流與表面夾角θc是來流方向與x 軸的夾角[24]。因此,當地壓力系數和剪切力系數分別為

各符號取值如下:θc滿足cos θc>0;σn、στ約為0.8;S 約為3~14,180~200 km 高度上,可取S=8。

采用當地化方法得到氣動力的壓力系數和剪切力系數公式,對應氣動力和氣動力矩為

當地微元總的氣動力為其相加之和:

2 氣動舵機輔助控制的超低軌衛星構型設計

根據整星初步構型和主要設備的分布方式,建立衛星模型如圖1 所示。主體部分橫截面為平行四邊形,且在其正x 端面中心位置有中心空腔,為設想采用的吸氣式電推進進氣道,同時具有降低前向氣動阻力的作用。在星上載荷與整星平臺間采用隔震連接,從而降低超低軌氣動環境及推進系統引入的振動影響。共4 幅氣動舵機分布于整星結構四邊位置,每塊舵機翼板由主翼和副翼組成,主翼轉軸過整星質心,副翼轉軸設置于主翼上,衛星質心位于其幾何中心。整星初步設計尺寸如圖2所示。

圖1 超低軌衛星氣動構型示意圖Fig.1 Pneumatic structure of ultra-low-orbit satellite

圖2 中,垂直舵機在星體坐標系xz 平面內沿x 軸對稱分布,其型心連線過衛星質心且相對質心對稱分布,主翼轉動角度α 與副翼轉動角度β 的定義如圖3 所示。

圖2 超低軌衛星初步設計尺寸Fig.2 Preliminary design dimension of ultra-low-orbit satellite

圖3 超低軌衛星主副翼轉角定義示意圖Fig.3 Definition of rotation angles of main wings and ailerons of ultra-low-orbit satellite

根據圖3 構型可知,當主翼帶動副翼同時轉動時(即β=0),僅產生過衛星質心的推力,而不會造成星體旋轉;與之相對,當副翼存在轉動角度時(即β≠0),翼面引入的氣動合力將偏離質心位置,從而形成轉動力矩,繼而改變衛星姿態。

3 氣動特性的仿真結果及分析

氣動力和力矩的計算結果受到很多因素的影響,結合衛星在軌運行過程中實際狀況及姿軌控制的需要,針對以下幾種情況進行計算分析:①在姿態控制過程中,衛星的姿態會有較大變化,因此計算氣動力和力矩隨衛星姿態角的變化情況;②利用氣動力和力矩輔助衛星的姿軌控制是通過氣動翼的偏轉實現的,因此計算氣動力和力矩隨氣動翼轉角的變化情況;③衛星在軌運行過程中,來流方向不斷改變,導致衛星所受氣動力和力矩發生變化,因此計算氣動力和力矩隨來流方向的變化情況。其他因素對氣動力和力矩的影響較為單一且變化范圍較小,因此設置為定值。仿真設計過程中初始狀態參數設置如表1 所示。

表1 仿真參數設置Table 1 Simulation parameters setting

1)姿態角變化

與一般工程實際的應用方式一致,衛星本體坐標系相對軌道坐標系的坐標轉換矩陣,采用歐拉角進行描述,按照z-y-x(3-2-1)順序旋轉;繞xb軸轉動角度φ,定義為滾轉角;繞yb軸轉動角度θ,定義為俯仰角;繞zb軸轉動角度Ψ,定義為偏航角。

先計算在各氣動力轉角均為零、來流方向指向x 軸負方向的情況下,氣動力和力矩隨衛星俯仰角和偏航角的變化。由于在本體坐標系下,當來流向量沿x 軸負方向時,單獨滾轉角的變化對衛星氣動受力沒有影響,計算的只是氣動力和力矩隨偏航角及俯仰角的變化情況。而實際上,由于來流方向與x 軸存在夾角,滾轉角變化也會帶來氣動力的微弱變化,但衛星的構型沿x 軸對稱,而且來流方向與x 軸方向的夾角很小,滾轉角變化導致的氣動力變化較小。

仿真沿各軸向氣動力和氣動力矩隨姿態角度的變化情況,如圖4~圖6 所示。在氣動力方面,如圖4(a)所示,Fx恒小于零;如圖5(a)所示,Fy主要受偏航角影響,且與偏航角呈正相關,受俯仰角影響較?。蝗鐖D6(a)所示,Fz主要受俯仰角影響,當偏航角較小時,與俯仰角呈負相關。三軸的氣動力量級均為10-1N。氣動力矩方面,如圖4(b)所示,當姿態角較大時,Mx變化較大;如圖5(b)所示,My主要受俯仰角影響,符號與俯仰角一致;如圖6(b)所示,Mz符號與俯仰角一致,且當各角較大時變化較大。Mx量級為10-2N·m,My、Mz量級為10-1N·m。

圖4 x 軸氣動力及力矩隨姿態角變化Fig.4 Change of aerodynamic force and moment of x-axis with satellite attitude

圖5 y 軸氣動力及力矩隨姿態角變化Fig.5 Change of aerodynamic force and moment of y-axis with satellite attitude

圖6 z 軸氣動力及力矩隨姿態角變化Fig.6 Change of aerodynamic force and moment of z-axis with satellite attitude

為了對氣動力矩受姿態角影響進行更準確地分析,截取部分仿真切面,分別在偏航角和俯仰角為零時,變動另一姿態角,對受到的氣動力矩進行分析,結果如圖7 所示。

圖7 星體姿態角度變化對氣動力矩的影響Fig.7 Change of aerodynamic moment with satellite attitude

從分析結果可知,偏航角變化影響的主要是Mz,其他軸力矩在零附近變化。俯仰角變化影響的主要是My,其他軸力矩在零附近變化。上述力矩變化量級在10-3~10-2N·m。

2)來流方向變化

為簡化計算過程,仿真中以極地軌道為參考,此時來流向量x 軸和z 軸分量始終不變,y 軸分量會隨著衛星運動呈現周期性變化。來流向量的y軸分量在-0.09~0.09 范圍內變化,對應來流向量與軌道坐標系x 軸負方向的夾角范圍為-5.14°~5.14°。設置3 個姿態角為零,各氣動翼轉角為零,根據來流向量的變化計算氣動力和力矩的變化規律,結果如圖8 所示。

圖8 不同來流方向下氣動力與氣動力矩Fig.8 Change of aerodynamic force and moment with air flow direction

氣動力方面,Fx始終為負,來流偏轉越大,其值越大,量級為10-1N;Fy與來流y 軸分量大致呈線性正相關,量級為10-2N;Fz為零。氣動力矩方面,Mx為負且變化很小,My在零附近輕微浮動,Mz正負與來流向量y 軸分量相反,量級為10-3~10-2N·m。

3)氣動舵機偏轉

在實際運用到氣動力和力矩時,更多是各姿態角為零而氣動翼轉角不為零,因此設置3 個姿態角為零,令來流指向x 軸負方向,計算氣動力和力矩隨氣動翼轉角的變化情況。偏轉主翼面,計算氣動力隨主翼面轉角的變化規律,結果如圖9 所示;偏轉副翼面,計算氣動力矩隨副翼面轉角的變化規律,結果如圖10 所示。

圖9 氣動力隨主翼面轉角的變化Fig.9 Change of aerodynamic force with rotation angle of main wings

圖10 氣動力矩隨副翼面轉角的變化Fig.10 Change of aerodynamic moment with rotation angle of ailerons

當主翼偏轉時,氣動阻力Fx始終為負,且隨其絕對值偏轉角增大而增大,量級為10-1N;側向力量級為10-2N,水平主翼面轉角為正時Fz為正,轉角為負時Fz為負;垂直主翼面轉角為正時Fy為負,轉角為負時Fy為正,可以利用這一特點通過偏轉主翼產生所需氣動力以輔助相對軌道位置維持。

當副翼偏轉時,相應軸的氣動力矩隨之變化,可以達到10-1N·m 量級,而通常動量輪組只能產生10-2~10-1N·m 的力矩,氣動力矩可以對姿態運動產生很大的影響,其姿態調整能力與傳統動量輪調姿方式相當,因此完全可將其應用于姿態穩定控制中。水平副翼偏轉導致My變化,二者大致線性正相關;垂直副翼偏轉導致Mz變化,二者大致線性正相關,可以利用這一特點通過偏轉副翼產生所需氣動力矩以輔助姿態控制。

由上述分析可知,分別偏轉水平副翼和垂直副翼,產生y 軸和z 軸氣動力矩,由于氣動力矩較大,可以為這兩軸姿態控制提供控制力矩,相應動量輪則不工作,因此在此過程中y、z 軸動量輪角動量積累為零。由于x 軸氣動力矩是由側向力提供的,相較于另外兩軸氣動力矩有數量級的差距,本文不考慮利用x 軸氣動力矩輔助姿態控制。

4 氣動舵機輔助的姿態控制策略設計與仿真

4.1 氣動舵機輔助的姿態控制策略設計

由第3 節氣動特性仿真分析結果可知,當副翼偏轉時,相應軸的氣動力矩隨之變化,而且在一定偏轉范圍內二者大致線性正相關。可利用這一關系,通過副翼偏轉提供所需氣動力矩,輔助動量輪進行姿態穩定控制。

基于氣動舵機輔助的姿態控制框圖如圖11 所示。具體控制過程如下:

圖11 采用氣動舵機輔助的姿態控制框圖Fig.11 Block diagram of attitude control with pneumatic steering gear

1)通過PID 控制律,根據期望姿態與實際姿態的偏差求得整星所需總的控制力矩Tc_total,并通過控制分配得到需要副翼提供的氣動力矩Tca和動量輪提供的控制力矩Tc。其中,動量輪控制部分采用傳統方式,不再贅述。

2)利用所需的氣動力矩Tca求得副翼所需偏轉角度β=βaTca/Ta,實際控制過程中通過控制副翼偏轉角度進行氣動舵機輔助姿態控制。其中,βa為氣動力矩和副翼轉角大致線性相關的邊界角度,根據實際氣動力矩曲線確定,Ta為βa對應的氣動力矩。

3)令相應軸的副翼偏轉角度為β,代入到氣動計算函數中,可得到副翼偏轉角度為β 時衛星受到的實際氣動力矩。

假設衛星為對地三軸穩定航天器,姿態動力學可用如下方程描述:

式中:Tx、Ty、Tz為包括控制力矩及衛星受到的其他干擾力;ωo為衛星軌道角速度。

式中:Kp、Ki、Kd分別為比例系數、積分系數、微分系數;Δφ、Δθ、ΔΨ分別為當前姿態角與目標姿態角的差值;Jx、Jy、Jz分別為各軸慣量。

通過式(15)計算得到期望控制力矩Tc_total,之后x 軸通過動量輪組響應,y、z 軸通過舵面偏轉響應,產生各軸實際控制力矩。調整選取合適的控制參數,分別取Kp={5,2.8,5},Kd={140,148,210}, Ki={0.06,0.06,0.06}。

在充分考慮系統誤差和環境干擾的前提下,可實現衛星姿態高精度、高動態范圍的控制。結合傳統高分辨率中低軌遙感衛星姿態控制精度要求,本文所設計基于氣動舵機輔助的超低軌衛星滿足如下姿態控制精度:①控制指向精度:≤0.005°(三軸3σ);②姿態穩定度:≤0.000 5(°)/s(三軸3σ)。

下面對氣動力矩隨副翼面轉角的變化規律進行綜合分析,實現利用所需的氣動力矩Tca求得副翼所需偏轉角度β。

計算主翼偏轉不同角度時氣動力矩隨副翼面轉角的變化情況,分別針對主翼面轉角為正和為負作圖,如圖12 所示。

圖12 主翼偏轉不同角度時氣動力矩隨副翼面轉角的變化Fig.12 Change of aerodynamic moment with rotation angle of ailerons for different rotation angles of main wings

結合主翼面轉角為零時氣動力矩隨副翼面轉角的變化規律,得到表2。

表2 氣動力矩與副翼偏轉對照關系Table 2 Relationship between aerodynamic moment and rotation angle of ailerons

上述分析是針對垂直翼面進行計算得到的,由此前分析可知,水平和垂直翼面的偏轉導致的氣動變化規律基本一致,因此利用水平翼面偏轉輔助姿態穩定控制時各項參數選取與此一致。

4.2 氣動舵機輔助的姿態控制策略仿真結果及分析

開展仿真分析工作的仿真環境由SolidWorks、ANSYS、MATLAB 等軟件組成,主要仿真步驟如下:

步驟1 基于設計的衛星構型,通過SolidWorks軟件進行三維建模。

步驟2 將三維模型導入ANSYS 軟件進行有限元網格劃分,進行微元特征計算。

步驟3 將微元特征數據文件導入MATLAB軟件進行氣動力和氣動力矩計算分析,并依此開展基于氣動輔助的姿態控制策略設計。

步驟4 通過MATLAB 軟件建立衛星姿態控制數學模型,本節主要開展此部分工作,仿真工況設置為通過副翼偏轉產生的氣動力矩進行y 軸和z 軸姿態控制,通過動量輪進行x 軸姿態控制。

步驟5 通過對仿真結果進行分析,對所設計控制策略的有效性進行驗證。

開展氣動舵機輔助姿態控制策略仿真分析的主要流程如圖13 所示。

圖13 仿真流程Fig.13 Flowchart of simulation

衛星運動的數學模型包含大量的微分方程,因此總的思路為采用四階龍格庫塔法求解微分方程。

根據以往衛星設計經驗,初步設置仿真參數及初始值如表3 所示。

表3 仿真參數說明Table 3 Description of simulation parameters

動量輪控制周期及舵機控制周期均為0.2 s。y 軸和z 軸通過舵機偏轉產生氣動力矩進行姿態控制,設置舵機偏轉角度上限為34.5°,x 軸由動量輪進行姿態控制。

由于在衛星實際運行過程中,還存在各類測量誤差及控制偏差,為使仿真過程盡可能接近實際情況,在可能出現誤差的地方設置隨機干擾,如表4 所示。

表4 仿真加入的隨機干擾Table 4 Random disturbance in simulations

在上述情況下,通過對姿態角控制的仿真,獲得氣動舵機輔助下衛星姿態角和角速度變化曲線分別如圖14 和圖15 所示。

圖14 氣動舵機輔助下的衛星姿態穩定過程Fig.14 Stabilization process of satellite attitude under assistance of pneumatic steering gear

圖15 氣動舵機輔助下的衛星姿態角速度穩定過程Fig.15 Stabilization process of satellite angular velocity under assistance of pneumatic steering gear

計算得到三軸的控制指向精度為[0.000 7°,0.003 3°,0.001 9°];姿態穩定度為[0.000 06,0.000 66,0.000 43](°)/s,與預想的姿態控制要求相差不大。

氣動舵機偏轉角度仿真結果如圖16 所示,舵機氣動力矩、氣動干擾力矩和動量輪輸出力矩的仿真結果如圖17~圖19 所示。

圖16 氣動舵機偏轉角度Fig.16 Deflection angle of pneumatic steering gear

圖17 舵機氣動力矩Fig.17 Aerodynamic torque of steering gear

圖18 舵機氣動干擾力矩Fig.18 Aerodynamic disturbance torque of steering gear

圖19 動量輪輸出力矩Fig.19 Output torque of momentum wheel

動量輪輸出力矩在y、z 軸分量均為0,只有x軸分量不為0,舵機偏轉分別在y、z 軸產生氣動控制力矩,上述結果表明成功實現x 軸動量輪組控制,y、z 軸氣動控制的功能。

5 結 論

1)本文研究了基于氣動舵機輔助的超低軌衛星姿態控制策略,并進行了計算仿真。利用稀薄氣體動力學中的自由分子流理論對氣動力與氣動力矩進行了建模;結合超低軌大氣環境特點和氣動舵機輔助需求,設計構建了超低軌氣動衛星構型;在此基礎上,針對不同姿態角、來流方向和氣動翼偏轉角度下的氣動力和氣動力矩進行了數值計算和仿真分析,確定了在當前設計下氣動力可達10-1N量級、氣動力矩可達10-1N·m 量級,其姿態調整能力與傳統動量輪調姿方式相當,充分驗證了氣動舵機輔助超低軌衛星姿態控制策略的可行性。

2)開展了基于氣動力輔助的姿態控制策略設計及仿真分析,數值結果表明,超低軌衛星的三軸指向精度可達[0.000 7°,0.003 3°,0.001 9°],姿態穩定度 可 達[0.000 06,0.000 66,0.000 43](°)/s。上 述 結 果已滿足傳統高分辨率中低軌遙感衛星的姿態控制精度要求,同時也驗證了基于氣動舵機輔助的超低軌衛星姿態控制策略具有差異化的技術先進性和廣闊的應用前景。

在未來研究中,將基于已有研究成果,深化論證方案的工程可實現方法,同時進一步探索基于氣動舵機輔助的超低軌衛星軌道機動方法,將超低軌衛星姿軌控全流程與氣動舵機輔助緊耦合,從而為實現高分辨新體制遙感衛星提供有力的技術支撐。

主站蜘蛛池模板: 国产高清毛片| 日韩av在线直播| 国产精品开放后亚洲| 精品国产自| 国产成人精品一区二区免费看京| 手机精品福利在线观看| 欧美一区精品| 女同久久精品国产99国| 一区二区午夜| 秋霞一区二区三区| 999福利激情视频| 在线免费a视频| 强乱中文字幕在线播放不卡| 成·人免费午夜无码视频在线观看| 亚洲男人的天堂在线| 99在线视频网站| 国产手机在线小视频免费观看| 国内精品视频区在线2021| 亚洲成人在线免费观看| 国产高清在线精品一区二区三区| 高清无码手机在线观看| 欧美成人在线免费| 日韩午夜片| 亚洲色图欧美激情| 日韩专区欧美| 四虎精品黑人视频| 国产v精品成人免费视频71pao| 国内丰满少妇猛烈精品播| 蜜芽一区二区国产精品| 99re这里只有国产中文精品国产精品| 国模私拍一区二区| 成人福利在线观看| www.99在线观看| 91蜜芽尤物福利在线观看| 国产在线第二页| a色毛片免费视频| 久久影院一区二区h| 亚洲黄色网站视频| 国产精品无码作爱| 一级爆乳无码av| 国产福利一区二区在线观看| 天天躁日日躁狠狠躁中文字幕| 夜色爽爽影院18禁妓女影院| 91精品国产综合久久不国产大片| 亚洲欧美精品一中文字幕| 国产三级毛片| 国产18在线| 色老头综合网| 青青青视频91在线 | 中文字幕亚洲精品2页| 青青草欧美| 欧美a在线看| 欧美亚洲激情| 亚洲经典在线中文字幕| www亚洲天堂| 毛片在线看网站| 欧美在线一二区| 欧美中文字幕一区| 国产91丝袜在线播放动漫| 日韩人妻无码制服丝袜视频| 波多野结衣的av一区二区三区| 国产精品播放| 2020精品极品国产色在线观看 | 日韩区欧美国产区在线观看| 精品国产女同疯狂摩擦2| 在线精品欧美日韩| 国产剧情国内精品原创| 日本手机在线视频| 亚洲三级成人| 久久国产黑丝袜视频| 999国产精品| 国产精品一区在线麻豆| 99久久精品免费视频| 国产在线麻豆波多野结衣| 亚洲第七页| 亚洲欧美在线精品一区二区| 91精品国产麻豆国产自产在线| 日韩在线网址| 亚洲精品国产精品乱码不卞| 精品国产成人a在线观看| 9久久伊人精品综合| 午夜精品一区二区蜜桃|