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

單太陽翼放氣干擾力矩分析

2018-02-28 00:43:33邊志強王鑫栗雙嶺董瑤海沈毅力曾擎洪振強宋效正
航天器工程 2018年1期

邊志強 王鑫 栗雙嶺 董瑤海 沈毅力 曾擎 洪振強 宋效正

(1 上海衛星工程研究所,上海 201109)(2 中國衛星海上測控部,江蘇江陰 214431)(3 上海航天技術研究院,上海 201109)

放氣現象在衛星中普遍存在,如微波部件為防止空間低氣壓放電一般都設計放氣孔,衛星推進管路在發射過程中會根據程序控制打開閥門進行管路放氣[1],衛星結構中廣泛應用的蜂窩夾層板、隔熱材料、膠黏劑等材料在真空、高低溫、粒子輻射等環境中蒸發、升華和分解釋放出氣體,這不僅影響材料本身的性能,還會對光學儀器等造成污染,甚至還會使姿態產生顯著變化。

目前,國內外都有針對衛星在軌干擾力矩辨識的研究,文獻[2-3]中以“國際空間站”為例,分析了飛船壓力艙壁受微小隕石和空間碎片撞擊出現漏氣小孔,提出漏氣小孔的噴氣產生干擾力矩,引起飛船姿態顯著變化,利用非線性濾波方法估計干擾力矩大小,但是該方法需要精確知道干擾力矩初值,而實際在軌時初值是很難精確得到的。文獻[4]針對衛星穩態下的運動,提出了一種基于特征系統實現算法的衛星周期干擾力矩估計方法,但是該方法僅應用于衛星穩態飛行時干擾力矩的估計;文獻[5]中針對軸對稱衛星提出了一種利用發動機噴氣過程角速度變化解算干擾力矩的方法,該方法在角速度較小時效果比較好,而且僅針對軸對稱航天器。

太陽翼的放氣模型很難建模,并且地面無法測試擾動力矩大小。衛星在軌運行時,太陽翼放氣干擾力矩很難直接測量,可通過測量其他物理參數進行間接估計。在軌實際估計太陽翼放氣干擾力矩,一方面可以及時修正控制參數,避免姿態頻繁振蕩,減少推力器頻繁工作和燃料過多消耗;另一方面,在軌估計數據可以作為后續衛星姿軌控設計輸入,進行推力器控制設計時可以放寬極限環閾值,根據放氣程度逐漸減弱采用閾值更小的控制極限環。本文首先給出太陽翼放氣以及影響衛星姿態的機理,推導利用衛星姿態動力學模型、姿態變化和推力器等參數估計太陽翼放氣產生的干擾力矩過程,最后以風云四號(FY-4)衛星轉移軌道段太陽翼對日定向過程為例,進行了理論模型的驗證,在軌實時估計太陽翼放氣力矩。

1 太陽翼在軌放氣機理

文獻[6-9]的相關研究表明,真空下氣體放氣是隨時間以負指數衰減,放氣率與時間關系可以描述為q=A·exp(-t/B),其中q為材料放氣率,A、B為放氣時間常數。文獻[10]中根據實驗數據擬合,得到的放氣速度表達式與理論分析的放氣負指數函數一致,說明了放氣速率過程的負指數函數模型的正確性。具體的放氣指數模型不是本文研究重點,不再進行深入研究。

太陽翼基板為碳纖維復合材料網格面板/鋁蜂窩芯夾層剛性結構,由碳纖維復合材料網格面板、鋁蜂窩芯、鉸鏈支座、聚酰亞胺薄膜等組成。蜂窩芯夾層由面板、蜂窩芯子和膠黏劑構成[11](見圖1)。由于地面制造和儲存使蜂窩夾層內含大量空氣,在低于0.01 Pa的空間環境下,蜂窩夾層內壓力增大,在材料表面上吸附的氣體從表面脫附,溶解于材料內部的氣體向真空邊界擴散,最后在界面上釋放、脫離材料,氣體分子不斷從其表面脫離釋放出來并揮發逸出,最終達到蜂窩內腔和外界的壓力平衡。根據理論模型可知,在放氣初階段放氣速率較快,但隨著時間衰減速率逐漸減小,趨于平穩。在發射主動段,星體、太陽翼也進行放氣,但主動段時外界有一定大氣壓力,對放氣速率有阻礙作用,且主動段不具備放氣速率和干擾力矩分析的條件。本文只研究星箭分離后,太陽翼對日定向過程中的放氣影響。

圖1 太陽翼蜂窩芯夾層結構示意圖Fig.1 Honeycomb sandwich interface schematic of solar wing

2 太陽翼放氣干擾力矩估計

估計靜止軌道衛星轉移軌道段太陽翼放氣干擾力矩時,空間環境干擾力矩很多,如重力梯度力矩、太陽光壓力矩、推力器羽流等,這些對太陽翼放氣干擾力矩的估計影響較小,分析如下:

(1)重力梯度力矩是周期性的,而本文估計的在軌放氣干擾力矩無周期性,此重力梯度力矩可以不考慮;推力器羽流影響只發生在推力器噴氣時刻,而在非噴氣時是對星體無干擾的,與在軌干擾力矩長時間持續存在現象不一致,可以忽略羽流影響。

(2)非對稱構型衛星的單太陽翼放氣力矩對衛星姿態變化影響較大,而雙太陽翼的對稱結構,放氣力矩為力偶模式,部分相互抵消,不容易被測量估計。本文只針對單太陽翼衛星的放氣情況進行分析。

(3)衛星在軌飛行過程中,姿態受到太陽光壓力矩、重力梯度力矩和氣動力矩的影響,橢圓軌道和靜止軌道衛星軌道特性決定了氣動力矩、重力梯度力矩方向不斷隨軌道周期性變化,并且持續時間較短。因此,在估計太陽翼放氣力矩時,可忽略重力梯度力矩和氣動力矩的影響。

作用在整個航天器上的光壓力矩Msun[12]為

Msun=-?Sρsuncosθ[(1-η)L+2ηcosθn]dS

(1)

式中:L為輻射源方向的單位矢量;0≤η≤1為航天器表面的反射系數;n為受照表面的法向單位矢量;ρsun是受照面處的光壓強度,ρsun=I0/c(I0為太陽輻射通量,c為光速);θ是受照面法向單位矢量與輻射源方向單位矢量之間的夾角;積分區域S表示遍及航天器承受光壓力的表面部分。

以靜止軌道衛星為例,取I0=1358 J/m2,η=0.9,θ=0°,積分有效區域S取(17~30)m2。則可計算出太陽翼對日定向過程中,單太陽翼的光壓力矩Msun值在0.2 mN·m~0.4 mN·m范圍,這個量級的干擾對衛星姿態影響較小,在進行放氣干擾力矩估計時可以忽略太陽光壓力矩的影響。

本文的衛星本體坐標系Ob-XbYbZb定義為:坐標原點Ob(衛星質心),ObXb軸、ObYb軸和ObZb軸為衛星的3個幾何軸;太陽翼本體坐標系Op-XpYpZp定義為:坐標原點Op(太陽翼安裝點),OpXp軸、OpYp軸和OpZp軸為太陽翼本體的3個幾何軸(見圖2)。衛星在對日定向巡航階段時,太陽翼不轉動,太陽翼電池片貼片面的法線與太陽光入射方向恰好相反。太陽翼基板的放氣孔在電池片貼片面的背面,放氣產生的反作用力與太陽翼電池片貼片面的法線一致,結合整星質心位置,可知太陽翼放氣作用力產生的擾動力矩在衛星本體坐標系+Xb軸方向,對衛星+Xb軸姿態角速度和姿態角有較大影響。

圖2 衛星太陽翼在軌放氣方向示意圖Fig.2 Air-bleed of solar wing schematic

2.1 帶干擾力矩的衛星動力學方程

衛星在轉移軌道段對日巡航狀態下,太陽翼放氣產生的干擾力矩是長時間持續存在的,對姿態的擾動通過陀螺測量得到角速度實時數據,進而輸入到姿態軌道控制計算機中,解算姿態控制推力器的工作狀態從而進行控制。衛星姿態動力學方程為

(2)

式中:J為衛星轉動慣量矩陣,ω為衛星慣性角速度,Mc為推力器控制力矩,Md為衛星受到的干擾力矩。姿態動力學方程中,轉動慣量可以在地面精確測量得到,姿控推力器控制力矩根據推力器安裝位置、推力大小和方向可以計算得到,衛星實時角速度利用陀螺可以測量得到。

(1)衛星轉動慣量:衛星在不同軌道階段、不同部件展開狀態下的轉動慣量,可在地面通過理論計算和試驗測試得到精準的數據,誤差在0.1%以內[13]。星箭分離后,用于姿態控制的燃料消耗較少,忽略整星質心的變化,認為此過程轉動慣量保持不變。

(2)推力器控制力矩:衛星巡航過程中的控制力矩由固定推力大小和方向不變的姿控推力器產生,Mc=L×F,L表示推力力臂,為衛星質心到推力作用點的向量,F為推力器工作推力。力臂L可以在地面精確地測量得到,且估計過程中保持不變。推力F在初始變軌階段認為保持不變。因此,控制力矩Mc可以較為準確的計算得到。

(3)衛星姿態角速度解算:衛星慣性角速度通過陀螺表頭構型以及各個接入系統陀螺的角速度測量輸出。根據陀螺表頭實際安裝指向,用矩陣A3×N描述陀螺安裝矩陣。實際在軌運行時,選取其中5個陀螺表頭的輸出值解算衛星慣性角速度,即

(3)

式中:A3×5為選取的5個表頭的安裝矩陣,ωxb、ωyb、ωzb為慣性系相對于衛星本體坐標系的三軸角速度,ωgi(i=1,2,3,4,5)為選取表頭的測量值,為慣性系相對于陀螺測量軸的三軸角速度。

2.2 太陽翼放氣干擾力矩計算

衛星姿態動力學方程表示如下:

(4)

式中:Jxx、Jyy、Jzz為衛星轉動慣量的主慣量,Jxy、Jzx、Jyz為衛星轉動慣量的慣性積;ωx、ωy、ωz為衛星三軸慣性角速度在本體系中的分量;Mcx、Mcy、Mcz為衛星在本體系下受到的三軸控制力矩;Mdx、Mdy、Mdz為衛星在本體系下受到的三軸干擾力矩。

慣量積相對主軸慣量較小,且經仿真計算知其影響很小可忽略不計,而角速度耦合項干擾力矩估計中不能忽略不計。因此,動力學方程簡化表示為

(5)

對式(4)[t0,t0+T]時間內進行積分(T為積分時間),對于其它坐標軸與本軸的耦合項,采用數值積分的方法,以陀螺角速度更新周期Tg(一般為0.5 s,滿足T=kTg,k為整數)為步長,在整個積分時長內進行積分求和。

(6)

控制力矩分為正向控制力矩和負向控制力矩,近似為恒定值。滾動軸控制力矩為Mcx+、Mcx-,推力器工作時長為Δtx+、Δtx-;俯仰軸控制力矩為Mcy+、Mcy-,推力器工作時長為Δty+、Δty-;偏航軸控制力矩為Mcz+、Mcz-,推力器工作時長為Δtz+、Δtz-。衛星遙測數據下傳各方向推力器工作時長為Δti、Δti-(i=x,y,z)。太陽翼放氣干擾力矩在整個過程中一直存在,而且變化緩慢。為了保證太陽翼放氣干擾力矩估計時效性,積分步長選取為秒量級的,認為在積分時長T內是不變化的。因此,在[t0,t0+T]時間內,有

(7)

因此,由式(5)~式(7)可以表示為

(8)

利用式(8)可以計算得到太陽翼放氣干擾力矩。

根據在軌數據需求確定得到迭代步長T1,每隔T1時長,取一個點作為本次積分的t0時刻,利用上述方法可以計算得到連續放氣干擾力矩。干擾力矩實時連續估計實施流程分為4步(見圖3):第一步,采集衛星角速度和推力器時間數據;第二步,確定相關時間參數;第三步,積分求和;第四步,解算放氣干擾力矩。

圖3 干擾力矩估計流程圖Fig.3 Flow chart of disturbance torque estimation

3 太陽翼放氣干擾力矩估計

3.1 理論仿真計算

以FY-4衛星參數為例,模擬轉移軌道段單太陽翼對日定向過程,估計太陽翼放氣干擾力矩。太陽翼展開后衛星慣量矩陣(單位:kg·m2)為

(9)

推力器產生Xb軸控制力矩Mcx+=Mcx-=18.9 N·m,產生Yb軸控制力矩Mcy+=Mcy-=19.75 N·m,產生Zb軸控制力矩Mcz+=Mcz-=15.4 N·m。陀螺角速度更新周期Tg為0.5 s,推力器工作計時更新周期以及積分步長T=16 s。為保證數據時效和數據量適中,迭代步長T1=16 s。衛星動力學模型中,光壓力矩模型采用式(1)模型和相關參數,并人為模擬了太陽翼放氣干擾力矩,其模型為Mdx=0.05·exp(-0.000 2·t)。仿真過程中采用與星上一致的開關線推力器控制規律,姿態角與姿態角速度的控制閾值分別為0.5°、0.004 5(°)/s。陀螺測量誤差為0.005(°)/s,太陽敏感器測量噪聲為0.1°。根據前述算法進行太陽翼放氣干擾力矩估計,估計結果如圖4所示。

圖4 放氣干擾力矩估計仿真Fig.4 Simulation result of disturbance torque

由圖4可知,太陽翼放氣干擾力矩按照負指數函數衰減,估計值與理論值誤差小于10%。仿真結果表明:本文算法能夠很好實現太陽翼放氣干擾力矩的實時估計,且太陽光壓力矩對太陽翼放氣干擾力矩估計基本無影響,在軌估計算法模型中不用考慮。

3.2 在軌實時估計

利用FY-4衛星星箭分離后在轉移軌道段的在軌數據,衛星參數按照3.1節選取。運用陀螺1、2、3、4、5的角度增量解算得到2016年12月10日三軸角速度如圖5~圖10所示(注:對日定向過程中,偏航軸不進行角度控制)。

圖5 慣性角速度ωx曲線Fig.5 Curve of inertial angular velocity ωx

圖6 慣性角速度ωx曲線(局部放大圖)Fig.6 Curve of inertial angular velocity ωx(amplify curve)

圖7 慣性角速度ωy曲線Fig.7 Curve of inertial angular velocity ωy

圖8 慣性角速度ωy曲線(局部放大圖)Fig.8 Curve of inertial angular velocity ωy(amplify curve

圖9 慣性角速度ωz曲線Fig.9 Curve of inertial angular velocity ωz

圖10 慣性角速度ωz曲線(局部放大圖)Fig.10 Curve of inertial angular velocity ωz(amplify curve)

2016年12月10日FY-4衛星太陽翼展開對日定向過程中,太陽翼放氣引起的衛星滾動、俯仰姿態角變化曲線如圖11所示,Xb軸方向姿態控制推力器工作計時曲線如圖12所示。結合圖11、圖12可知,在衛星對日定向過程中,滾動角曲線呈現單邊極限環,說明在Xb軸方向存在一個干擾力矩,即太陽翼放氣產生的+Xb軸方向干擾力矩,則衛星推力器產生-Xb軸方向的控制力矩。因此,產生-Xb軸控制力矩的推力器一直工作,計時不斷增加,而+Xb軸控制力矩的推力器不工作,推力器計時不變化。

圖11 太陽翼放氣引起的衛星滾動、俯仰姿態角變化曲線Fig.11 Change curve of attitude angles

圖12 ±Xb軸方向推力器計時曲線Fig.12 Jet timing data of thrusters in ±Xb axis direction

根據上述數據以及本文提出的方法,得到太陽翼放氣干擾力矩估計結果如圖13所示。

圖13 風云四號衛星太陽翼展開后放氣干擾力矩估計值、擬合函數及誤差Fig.13 Estimation result of disturbance torque caused, fitting function and fitting error of FY-4

將放氣干擾力矩按照指數函數進行曲線擬合,得到分段時間函數關系,并與在軌估計值作比較,得到擬合函數的誤差值。放氣干擾力矩分段時間函數為

(10)

根據在軌實時估計結果,可知:

(1)衛星慣性角速度ωx變化較快,推力器頻繁工作進行姿態控制,是因為受到太陽翼放氣干擾力矩的影響。在軌放氣干擾力矩大小,前0.64 h放氣干擾力矩按e-1.8t指數函數衰減,0.64 h后按e-0.485t指數函數衰減。這是因為隨著太陽翼蜂窩夾層內氣體不斷排出,內部壓力逐漸變小,氣體排出產生的干擾力矩衰減很快。

(2)根據單太陽翼構型可知,放氣干擾力矩作用在+Xb軸方向上,根據在軌陀螺數據和推力器工作時間計算太陽翼放氣產生的干擾力矩最大達到37 mN·m,比式(1)計算的單太陽翼太陽光壓力矩約大2個數量級,并且太陽翼放氣時間持續近11 h。

4 結束語

本文根據衛星姿態動力學原理,提出了太陽翼放氣干擾力矩估計方法,通過理論仿真,分析了方法的正確性,并運用FY-4衛星在軌實時估計驗證了放氣干擾力矩按照負指數模型衰減。該方法面向工程實際,不需要建立復雜的復合材料放氣模型,即可精確建立單太陽翼衛星姿態動力學特性。

衛星在軌運行受到太陽光壓力矩、重力梯度力矩和氣動力矩等干擾,且多具有周期性,能夠精確建模,可采取前饋控制或反饋控制對干擾進行吸收;而衛星轉移軌道段太陽翼放氣干擾力矩對姿態擾動較大,估計太陽翼放氣干擾力矩對衛星穩態控制有重要意義。

(1)利用估計值及時修正控制參數,避免姿態頻繁振蕩,減少推力器頻繁工作消耗過多燃料。

(2)可作為后續衛星姿軌控設計參考,考慮放氣干擾力矩影響,提高系統可靠性和安全性。

(3)在太陽翼制造過程中,適當增加或加大放氣孔,使在主動段或轉移軌道段初期放氣盡快結束,以減少對姿態的影響。

References)

[1] 王獻忠,劉赟,張麗敏,等.星載計算機及關鍵流程可靠性設計[C]//第四屆中國指揮控制大會論文集.北京: 中國指揮與控制學會,2016:492-496

Wang Xianzhong,Liu Yun,Zhang Limin,et al.Reliable design for on-board computer and key process[C]//The Fourth Conductor in China Control the General Assembly Talk About Collection of Essays.Beijing:The Fourth Conductor in China Control the General Assembly Talk About Collection of Essays.2016:492-496 (in Chinese)

[2] J W Kim,J L Crassidis,S R Vadali,et al.ISS leak localization using attitude response[J]. AIAA Guidance, Navigation, and Control Conference and Exhibit. Washington D.C.:AIAA,2001

[3] 李海軍,黃顯林,胡海東,等.基于濾波算法和飛船姿態響應尋找漏氣小孔[C]//第25屆中國控制會議論文集(上冊).哈爾濱:哈爾濱工業大學,2006:492-496

Li Haijun, Huang Xianlin, Hu Haidong, et al. Spaceship leak localization based on the attitude response and filtering algorithm[C]//The 25th Chinese Control Conference(up).Harbin:Harbin Institute of Technology,2006:492-496 (in Chinese)

[4] 雷靜,劉瑩瑩,周鳳岐,等.衛星在軌周期干擾力矩辨識與補償方法的研究[J].西安:西北工業大學學報,2009,27(3):396-400

Lei Jing,Liu Yingying,Zhou Fengqi,et al.A new method for identification of and compensation for periodically disturbing torques of on-orbit satellite[J].Xi’an:Journal of Northwestern Polytechnical University,2009,27(3):396-400 (in Chinese)

[5] 梁彤,張奕群.一種空間飛行器軌控發動機干擾力矩的測試方法[J].現代防御技術,2008,36(1):62-65

Liang Tong,Zhang Yiqun.Method of testing disturbing torque of divert thrusters for the spacecraft[J].Modern Defence Technology,2008,36(1):62-65 (in Chinese)

[6] 楊春光,肖尤明,陳楠,等.真空下非金屬材料放氣模型與研究綜述[J].真空,2016,12(3):48-53

Yang Chunguang,Xiao Youming,Chen Nan,et al.A comprehensive survey of development of outgassing model for nonmetal materials in vacuum[J].Vacuum,2016,12(3):48-53 (in Chinese)

[7] 羅艷,王魁波,張羅莎.聚合物的放氣分率與放氣模型研究[J].真空科學與技術學報.2015,9(9):1100-1115

Luo Yan,Wang Kuibo,Zhang Luosha.Modelling and characterization of polymer outgassing behavior[J].Chinese Journal of Vacuum Science and Technology,2015,9(9):1100-1115 (in Chinese)

[8] 張亞平,朱穎峰,劉湘云,等.基于擴散放氣模型的杜瓦真空壽命分析[J].真空科學與技術學報,2014,34(1):28-31

Zhang Yaping,Zhu Yingfeng,Liu Xiangyun,et al.Influence of diffusion out-gassing on life-time of dewar vessels[J].Chinese Journal of Vacuum Science and Technology,2014,34(1):28-31 (in Chinese)

[9] Schindler N,Schleurner D,Chr Edelmann.Measurement of partial outgassing Rates[J].Vacuum,1996,47(4):351-355

[10] 陳濤,李玉忠,許忠旭,等.真空熱試驗中材料放氣的放氣量及其導熱問題[J].航天器環境工程,2006,2(23):103-106

Chen Tao,Li Yuzhong,Xu Zhongxu et al.Material outgassing amount and thermal conduction in vacuum thermal test[J].Spacecraft Environment Engineering,2006,2(23):103-106 (in Chinese)

[11] 袁家軍.衛星結構與設計(下)[M].北京:宇航出版社,2004:216-218

Yuan Jiajun.Satellite structure and design(next)[M].Beijing:Astronautics Press,2004:216-218 (in Chinese)

[12] 劉善武,萬松,容建剛.航天器空間環境干擾力矩分析與仿真研究[J].航天控制,2015,2(33):78-90

Liu Shanwu,Wan Song,Rong Jiangang.The analysis and simulation of aircraft space environment disturbance torque[J].Aerospace Control,2015,2(33):78-90 (in Chinese)

[13] 徐小方,張華.飛行器轉動慣量測量方法研究[J].科學技術與工程,2009,3(6):1653-1660

Xu Xiaofang,Zhang Hua.Method of testing the aerocraft’s moment of inertia[J].Science Technology and Engineering,2009,3(6):1653-1660 (in Chinese)

主站蜘蛛池模板: 乱码国产乱码精品精在线播放| 国产在线91在线电影| 一级毛片中文字幕| 乱人伦99久久| 国产特一级毛片| 伊在人亚洲香蕉精品播放| 国产精品久久久久婷婷五月| 国产av剧情无码精品色午夜| 一本色道久久88亚洲综合| 全色黄大色大片免费久久老太| 国产爽歪歪免费视频在线观看| 亚洲浓毛av| 永久在线播放| 亚洲天堂精品视频| 日本午夜精品一本在线观看 | 99精品免费在线| 国产精品密蕾丝视频| 麻豆国产精品| 国产尤物jk自慰制服喷水| 欧美一区国产| 久久亚洲国产最新网站| AV网站中文| 亚洲综合狠狠| 国产精品区视频中文字幕| 久久伊人操| 日韩精品亚洲一区中文字幕| 精品久久久久久中文字幕女| 日韩性网站| 午夜天堂视频| 91娇喘视频| 精品成人一区二区三区电影| 91久久性奴调教国产免费| 欧美色伊人| 国产精品乱偷免费视频| 久青草免费在线视频| 亚洲二区视频| 亚洲成a人片7777| 午夜精品久久久久久久无码软件 | 亚洲国产精品VA在线看黑人| 国产在线观看成人91| 老汉色老汉首页a亚洲| 91精品国产自产在线老师啪l| 亚洲妓女综合网995久久| 国产成人久视频免费 | 国产靠逼视频| 婷五月综合| 国产一区亚洲一区| 久久天天躁夜夜躁狠狠| 国产一区亚洲一区| 国产内射在线观看| 色综合狠狠操| 国产不卡国语在线| 国产一区免费在线观看| 丁香五月激情图片| 国产人在线成免费视频| 免费人欧美成又黄又爽的视频| 四虎国产成人免费观看| 四虎成人精品| 国产在线一二三区| 欧美视频在线播放观看免费福利资源| 自偷自拍三级全三级视频| 国内嫩模私拍精品视频| 精品一區二區久久久久久久網站| 无码高潮喷水在线观看| 国产精品久线在线观看| 天天综合网色| 国产女人在线| 欧美日韩福利| 久久频这里精品99香蕉久网址| 狠狠色狠狠色综合久久第一次| 伊人网址在线| 精品剧情v国产在线观看| 久久永久免费人妻精品| 久久99蜜桃精品久久久久小说| 国产欧美另类| 精品久久香蕉国产线看观看gif| 欧美色图久久| 国产69精品久久久久孕妇大杂乱 | 国产二级毛片| 精品国产自在现线看久久| 国产精品无码AⅤ在线观看播放| 欧美亚洲日韩不卡在线在线观看|