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

航天器在軌全周期熱變形分析方法

2016-12-29 11:11:34劉國青羅文波童葉龍范立佳
航天器工程 2016年6期
關鍵詞:變形結構分析

劉國青 羅文波 童葉龍 范立佳

(北京空間飛行器總體設計部,北京 100094)

航天器在軌全周期熱變形分析方法

劉國青 羅文波 童葉龍 范立佳

(北京空間飛行器總體設計部,北京 100094)

提出一種適應于在軌全周期熱變形的分析方法,采用基于熱傳導算法進行“熱分析模型-結構分析模型”全周期溫度場映射,利用數學擬合算法開展對各類結果數據的分析,通過相關程序實現全周期多工況溫度場映射、計算、數據分析的自動化。對某遙感衛星進行全周期熱變形分析,結果表明:全周期溫度場映射時間由天縮短至小時量級,溫度場映射精度可控制在1%以內,相對于以往基于極端工況的熱變形分析方法,可顯著地提升分析精度與驗證覆蓋性,獲得在軌熱變形量級、全周期變化規律。文章的研究結果可為航天器熱穩定設計提供參考。

航天器;在軌全周期;熱變形;穩定性

1 引言

隨著對地、對天觀測航天器指標要求的日益提高,高圖像定位精度成為高性能遙感航天器的典型特征。在軌結構變形直接影響相機、星敏感器、陀螺等關鍵部件自身空間指向及彼此間的幾何關系,甚至影響相機內部各鏡片間的空間位置關系,是決定圖像定位精度、相機成像質量的重要因素之一[1-3]。結構在軌熱變形在相機安裝處引起的位移一般為微米級,對于低分辨率觀測航天器,這些擾動可以忽略;但對分辨率優于1 m的航天器,則必須考慮熱變形擾動影響,因為這些影響可能直接決定了成像質量能否達到設計指標。一般情況下,微米級的結構變形可能導致角秒級的設備安裝面法向指向變化,進而出現米級的成像誤差。同時,與其他因素相比,航天器結構在軌熱變形具有一定的隨機性,很難通過后期在軌處理消除其影響。因此,在軌穩定性對高分辨率航天器的性能指標至關重要,在地面研制階段就應結合航天器系統需求開展航天器高穩定結構設計、驗證工作。

仿真分析是航天器機械系統研制的有力支撐,對于設計工況復雜、影響因素眾多的結構熱穩定設計而言,在軌熱變形分析直接支撐了系統指標分配、結構研制、熱變形地面試驗方案設計等各研制環節[4-5],國內已經開展的結構熱變形分析工作多數是針對結構熱膨脹、吸濕性的理論研究[6-8],以及零膨脹鋪層設計的研究[9]等,對航天器結構的在軌熱變形分析及試驗驗證主要立足于某特定溫度場或模擬溫度場,如模擬在軌工況的最高溫工況、最低溫工況、最大溫差工況等,或通過施加最大包絡載荷實現對在軌熱變形的預估及結構低膨脹設計[10]。隨著結構熱穩定性要求的日益提高,國外對結構熱穩定性的研究已經深入到在軌微裂紋、微蠕變等領域[11-13],并開展了結構熱穩定高保真仿真方法研究、影響因素靈敏度研究[14-16]等,這些研究的基礎在于對在軌全周期熱變形的高效、高精度仿真。在軌全周期熱變形分析還能為在軌成像標定策略的制定提供參考,進而實現在傳統圖像修正方法的基礎上進一步提升修正精度。因此,新一代遙感航天器研制對在軌全周期熱變形分析需求日益迫切。由于在軌全周期熱變形分析涉及的工況數量較大,一般可達數百、甚至上千個溫度工況,且面臨機熱耦合效應復雜、溫度邊界及力邊界模擬難度大等技術瓶頸,因此對熱控設計溫度場數據與結構分析溫度場數據高效、高精度映射,以及多工況下熱變形高保真分析、結果高效處理及判讀等,均提出了更高要求。在ESA發布的部分設計資料里,雖然提到了熱變形分析方法[17],但沒有闡述如何開展在軌全周期各個時刻的熱變形分析。目前,國內關于全周期熱變形分析的案例較少,支撐全周期熱變形分析、實現海量數據快速映射的方法更是鮮有提及;雖然有學者依據熱分析模型重新劃分結構分析模型,然后從數據文件中讀取相應節點的各個時刻的溫度數值,按照不同載荷工況的形式寫入計算文件進行溫度場分析[18],但是數據轉換過程較多,影響了全周期熱變形分析的效率。

本文首先對適應于在軌全周期熱變形分析方法進行探討,重點介紹了基于熱傳導算法的“熱分析模型-結構分析模型”全周期溫度場映射方法,以及基于數學擬合算法開展結果數據處理及判讀方法,闡述了全周期熱變形分析流程。基于上述方法對某遙感衛星進行全周期熱變形算例分析,獲取了在軌熱變形量級、全周期變化規律等,并與傳統分析方法進行了對比,結果表明本文提出的方法可顯著地提升分析精度與驗證覆蓋性。

2 全周期熱變形分析方法

2.1 總體思路

用于熱變形分析輸入的溫度場,通常是基于Thermal Desktop、I-DEAS、UG等軟件開展的熱分析得到,用于熱分析的數學模型與用于結構分析的數學模型(一般為通過PATRAN、ABAQUS等有限元分析軟件建立的結構分析有限元模型)在節點位置、網格離散程度、建模簡化方式等方面均存在差異性。例如,某航天器高穩定載荷適配支撐結構熱分析模型約有6000個節點,結構分析模型有16 000個節點[10]。因熱分析與結構分析所采用的軟件差異性,以及二者分析模型的差異性,在開展熱變形分析前,首先要將熱分析溫度場映射至結構分析模型上,進而實現溫度數據從熱分析模型傳遞至結構分析模型,并作為結構分析輸入載荷。本文提出的方法是首先生成無溫度場的結構分析計算文件,然后進行全周期工況判讀及分析,基于熱傳導算法實現熱分析模型向結構分析模型的溫度場映射,并對可能存在的不能映射節點和奇異節點(溫度遠高于或遠低于在軌實際溫度的節點)進行二次映射,繼而生成映射后的溫度場及結構分析計算文件,判斷無誤后進行計算,并對全周期熱變形分析結果進行數據擬合、生成報告。具體流程如圖1所示。

圖1 全周期熱變形分析流程Fig.1 Flow chart of whole cycle thermal deformation analysis

2.2 基于熱傳導算法的全周期溫度場映射

2.2.1 熱傳導算法概述

從映射算法上講,以往熱變形分析主要采用基于幾何差值的映射算法,該算法僅與空間位置相關,是目前商業軟件中廣為采用的映射方法。其局限性在于,對于非連續結構、不同組件連接結構的映射工況,易產生映射奇異的現象,無法識別各部位或各組件之間溫度的差異性。從溫度場映射實現工具來說,目前廣泛采用基于有限元商業軟件進行溫度場映射,此類方法對于溫度場單次映射較為通用,但對于全周期熱變形分析則具有一定的局限性,具體表現為:①溫度場導入及映射功能主要基于手動實現,很難滿足全周期成百上千個工況溫度場高效映射分析。②商業軟件一般僅內嵌基于幾何算法的映射方法。③因單位不同、設計師不同,熱分析過程可采用Thermal Desktop、I-DEAS、UG等不同軟件,進而導致熱分析結果數據格式存在顯著差異,商業軟件在數據導入、映射方式上對于各類熱分析軟件適應性較差。

本文提出的基于熱傳導算法的全周期溫度場映射方法,能以既有節點溫度場為基礎,依據結構熱傳導特性進行映射計算,避免因2個組件空間距離較近、但并不屬于同一溫度范圍的節點發生映射關系。同時,采用二次開發程序對映射方法進行封裝,分別建立適應于不同熱分析軟件的溫度場映射模塊。通過開發與有限元商業軟件前后處理工具的接口,實現溫度場映射批量處理、結果數據批量處理、結果數據批量判讀等,可顯著地提升映射精度和分析效率。通過圖2(a)一個簡單的結構組件分析模型,可以證明基于熱傳導算法進行溫度場映射的優勢。算例中的結構組件由高溫結構、低溫結構、室溫結構3部分組成,彼此間存在隔熱層,可以阻斷熱傳遞。基于幾何插值映射方法無法考慮低熱導率層的影響,見圖2(b);而基于熱傳導算法得到的溫度場映射結果,可更為真實地反映實際溫度情況,見圖2(c)。

圖2 基于幾何算法與基于熱傳導算法的溫度場映射結果對比
Fig.2 Comparison between mapping based on geometry and mapping based on thermal conductivity

2.2.2 溫度場映射過程

基于熱傳導算法的全周期溫度場映射過程包括3個步驟。

(1)構建熱分析模型單元節點和結構分析模型單元節點之間的對應關系,如圖3所示。

圖3 熱分析單元與結構分析單元溫度場映射對應關系Fig.3 Mapping relationship of temperature field between thermal analysis element and structural analysis element

(2)基于數學方法實現熱分析模型節點與結構分析模型節點的關聯,一般準則為:覆蓋熱分析模型節點的結構分析模型節點溫度,按照熱分析模型節點溫度取值。

使用有限元形函數獲取加權系數ai。

(1)

寫成矩陣形式為

Tt=ATf

(2)

式中:Tt為熱分析模型節點溫度矩陣;A為權重系數矩陣;Tf為結構分析模型節點溫度矩陣。

通過用熱控材料替代結構分析模型材料(例如用MAT4材料卡片替換MAT1材料卡片),基于結構分析模型可計算得到熱傳導矩陣Ct。

(3)求解如下的插值方程。

(3)

式中:q為拉格朗日乘子。

通過式(3)即可求解結構分析模型節點溫度矩陣Tf。

2.3 數據擬合算法

航天器在軌熱變形分析的目的,是獲取關鍵設備指向變化或各設備間的夾角變化,由此引申出采用何種方式來表征設備指向及其夾角的問題。目前,國內外廣泛采用的表征方式有2種。

(1)對于光學相機、星敏感器等設備,主鏡、次鏡等關鍵部件均沿設備軸向且近似在一條直線上,可選取此線上的多個關鍵點,應用“多點擬合線”的方式獲取設備指向。

(2)選取設備安裝面上的多個關鍵點,采用多點擬合面的形式獲取設備安裝面矢量,以此模擬設備安裝指向,此種方法對于各類設備均具有通用性。

由于結構分析結果一般為有限元模型節點位移,因此要借助其他程序并選取相應數學算法,對有限元分析結果進行二次處理和判讀[19-20]。

對于直線矢量計算,設待擬合直線矢量n個節點的坐標為(x1,y1,z1),(x2,y2,z2),…(xn,yn,zn),寫成如下矩陣形式。

(4)

計算式(4)的協方差矩陣D如下。

(5)

式中:

(6)

對于平面法線向量計算,設待擬合平面法線矢量的n個節點的坐標為(x1,y1,z1),(x2,y2,z2),…,(xn,yn,zn),寫成如下矩陣形式。

(7)

由式(7)各列減去各自的均值,得到矩陣R如下。

(8)

3 算例分析

3.1 溫度場映射

某遙感衛星3個星敏感器通過支架安裝于相機承力框上,相機安裝在衛星結構平臺上,其中整星機械坐標系Z向為衛星縱向(相機對地觀測方向),整星機械坐標系X向、Y向為衛星橫向(相機承力框及載荷適配結構面內方向)。該衛星每天運行15個軌道周期,熱控設計時要對衛星全生命周期所有極端工況取一個最大包絡,即衛星在軌運行每天承受的溫度工況均不會超過目前給定的15個軌道周期狀態,以載荷適配結構為例,其中3個典型位置的15個軌道周期熱分析節點溫度見圖4。

根據高定位精度設計需求,須開展15個軌道周期不同姿態下相機成像、數傳記錄等關鍵時刻點的熱變形分析,由此獲取相機安裝面法向轉角、星敏感器安裝面法向轉角,以及相機安裝面法向與星敏感器安裝面法向間的夾角變化。同時,分析不同的姿態、星敏感器和工作模式下上述各項分析結果的變化規律,為高定位精度指標分析提供支撐。

將采用Thermal Desktop軟件得到的全周期熱分析溫度場作為輸入,采用圖1仿真流程、基于熱傳導算法進行全周期近千余時間點溫度場映射,并對映射奇異節點進行二次修正,生成可用于NASTRAN軟件進行有限元分析的批處理求解文件。此外,通過MATLAB程序實現“熱分析溫度場輸入-溫度場映射-溫度場修正-有限元計算”高度集成化與自動化,進而實現全周期溫度場映射時間由天縮短至小時量級。

為驗證映射精確性,選取第一軌道周期溫度梯度較大4個時刻點,分別對應第一軌道周期溫度場“正弦曲線”的起點時刻、波谷時刻、波峰時刻、終點時刻。將映射前的熱分析溫度場與映射后的結構分析(有限元分析)溫度場進行對比,由圖5熱分析溫度場與結構分析溫度場對比結果可以看出,基于熱傳導算法可實現熱分析與結構分析溫度場的精確匹配,由熱分析模型到結構分析模型的溫度場映射精度可控制在1.00%以內(詳見表1)。

圖4 載荷適配結構3個典型位置的15個軌道周期熱分析節點溫度Fig.4 Node temperature of thermal analysis for three typical structural position through 15 orbit cycles

圖5 第一軌道周期典型時刻熱分析溫度場與結構分析溫度場對比
Fig.5 Comparison between temperature field of thermal analysis and structural analysis at typical time point of the first cycle

表1 溫度場映射精度

Table 1 Mapping precision of temperature field

時刻模型類型及誤差最高溫度點最低溫度點起點時刻熱分析模型38.49℃-51.78℃結構分析模型38.50℃-51.80℃誤差0.03%-0.04%波谷時刻熱分析模型38.22℃-31.37℃結構分析模型38.20℃-31.40℃誤差0.05%-0.10%波峰時刻熱分析模型69.90℃-34.30℃結構分析模型69.90℃-34.30℃誤差0.00%0.00%終點時刻熱分析模型38.61℃-51.16℃結構分析模型38.60℃-51.20℃誤差0.03%-0.08%

3.2 全周期熱變形分析結果

在以往熱變形評估工作中,主要選取模擬在軌工況的極端工況進行分析,尤其是以某個高溫或低溫狀態的均勻溫度作為輸入載荷,基于上述思路,針對圖5所示第一軌道周期典型時刻溫度場,選取4個極端溫度點(38.50 ℃,69.90 ℃,―34.30 ℃,―51.78 ℃)作為均勻溫度載荷輸入(模擬工況),開展熱變形分析,獲取星敏感器安裝面法向與相機安裝面法向夾角變化,并與全周期分析結果(真實工況)進行對比,見圖6。從圖6可以看出:4個模擬工況雖然可以反映一個周期內變形的平均值,但不能覆蓋全周期各時刻點可能出現的真實變形情況,說明采用全周期熱變形分析的必要性和優勢。

圖6 星敏感器安裝面法向與相機安裝面法向夾角變化對比Fig.6 Comparison of angle variation between plane normal line of one star sensor support and camera support

在軌運行期間,衛星有效載荷(如光學遙感衛星的對地或對天觀測相機)部分會采取精密控溫、相對常溫變化僅為幾攝氏度,而衛星平臺部分相對常溫存在幾十攝氏度的溫度波動。選取第一軌道周期變形量最大時刻點進行變形分析,由圖7~9整星、平臺、主承力立柱、相機適配支撐結構等部分變形云圖可以看出,在軌溫度交變引起的平臺部分變形在幾百微米,接近毫米級。而由圖10、圖11相機、相機主承力框變形云圖可以看出,相對于平臺部分而言,相機主承力框變形相對較小。

圖7 整星變形
Fig.7 Thermal deformation of satellite

圖8 主承力立柱變形
Fig.8 Thermal deformation of main truss

圖9 相機適配支撐結構變形
Fig.9 Thermal deformation of camera adaptor

圖10 相機部分變形
Fig.10 Thermal deformation of camera

圖11 相機承力框變形
Fig.11 Thermal deformation of camera support frame

通過全周期熱變形分析,獲取了1~15軌道周期3個星敏感器安裝面與相機安裝面法向夾角變化,見圖12。由分析結果可以看出:①3個星敏感器安裝面與相機安裝面法向夾角變化呈現正弦周期性變化,每一圈對應一個完整正弦波,此種變化趨勢與圖4所示的熱分析溫度場周期性變化情況相對應。②星敏感器安裝面與相機安裝面法向夾角變化存在顯著差異性,如+X+Y星敏感器安裝面與相機安裝面法向全周期最大夾角變化超過30″,而另外2個星敏感器安裝面與相機安裝面法向最大夾角變化未超過20″。③3個星敏感器支架同一變形形態出現時刻存在差異,即一個星敏感器支架變形位于“波峰”之時,另一個可能位于“波谷”。由此,獲得了不同星敏感器安裝面與相機安裝面法向最大夾角變化量級及彼此差異性,可為整星在軌穩定性評估、衛星定位精度評估等提供重要參考。

4 結論

本文提出了在軌全周期熱變形分析方法,基于熱傳導算法進行“熱分析模型-結構分析模型”溫度場映射,基于數學擬合算法開展結果數據處理,獲得了航天器在軌全周期熱變形量級、全周期變化規律等,可得出如下結論。

(1)基于熱傳導算法并通過相關程序實現全周期千余時間點映射過程的自動化,可將全周期溫度場映射時間由天縮短至小時,溫度場映射精度可控制在1%以內,相對于傳統分析方法顯著地提升了分析精度和全周期覆蓋性。

(2)獲取了整星在軌全周期熱變形情況,從變形云圖可以得到整星平臺、相機、星敏感器支架等關鍵部位變形狀態、量級及宏觀變形傳遞趨勢。

(3)通過各星敏感器安裝面與相機安裝面夾角變化全周期分析,可以看出3個星敏感器安裝面與相機安裝面夾角變化呈現正弦周期性變化,且變形量級、同一變形形態出現時刻均存在差異。

本文的熱變形分析方法及所獲取的變形規律具有一定的普適性,可用于航天器機械系統的熱穩定設計,對于高精度、高穩定性航天器的研制具有參考價值。

References)

[1]Uguen G,Luquet P,Chassat F. Design and development of the 2m resolution camera for Rocsat-2[C]//Proceedings of the 5th International Conference on Space Optics (ICSO 2004). Paris: ESA,2004: 173-180

[2]Lamard J L,Gaudin Delrieu C,Valentini D,et al. Design of the high resolution optical instrument for the Pleiades HR earth observation satellites [C]//Proceedings of the 5th International Conference on Space Optics (ICSO 2004).Paris: ESA,2004: 173-180

[3]Trigo J. Dimensional stability characterization of carbon fiber with epoxy and cyanate ester resin laminates due to moisture absorption [C]//Proceedings of Conference on Spacecraft Structures,Materials & Mechanical Testing. Paris: ESA,1996: 371-376

[4]譚維熾,胡金剛.航天器系統工程[M].北京:中國科學技術出版社,2009

Tan Weichi,Hu Jingang. Spacecraft systems engineering [M]. Beijing: China Science & Technology Press,2009 (in Chinese)

[5]彭成榮.航天器總體設計[M].北京:中國科學技術出版社,2011

Peng Chengrong. System design for spacecraft [M]. Beijing: China Science & Technology Press,2011 (in Chinese)

[6]袁家軍.衛星結構設計與分析[M].北京:中國宇航出版社,2004

Yuan Jiajun. Design and analysis of satellite structures [M]. Beijing: China Astronautics Press,2004 (in Chinese)

[7]Chen Liemin. Spacecraft structures and mechanisms [M]. Beijing: China Science and Technology Press,2009

[8]陳烈民,楊寶寧.復合材料的力學分析[M].北京:中國科學技術出版社,2008

Chen Liemin,Yang Baoning. Mechanical analysis for composite materials [M]. Beijing: China Science & Technology Press,2008 (in Chinese)

[9]武勇斌.微膨脹系數復合材料構件鋪層設計研究[D].哈爾濱:哈爾濱工業大學,2007

Wu Yongbin. Study and design on ply patterns of low thermal expansion of composites [D]. Harbin: Harbin Institute of Technology,2007 (in Chinese)

[10]劉國青,阮劍華,羅文波,等.高穩定結構熱變形分析與試驗驗證方法研究[J].航天器工程,2014,23(2):19-25

Liu Guoqing,Ruan Jianhua, Luo Wenbo,et al. Research on thermal deformation analysis and test verification method for spacecraft high-stability structure [J]. Spacecraft Engineering,2013,23(2): 19-25 (in Chinese)

[11]Smith P A. Comprehensive composite materials [M]. Amsterdam: Elsevier,2000

[12]Bathias C. An engineering point of view about fatigue of polymer matrix composite materials [J]. International Journal of Fatigue,2006,28: 1094-1099

[13]Pagano N J,Schoeppner G A,Kim R,et al. Steady-state cracking and edge effects in thermo-mechanical transverse cracking of cross-ply laminates [J]. Composites Science and Technology,1998,58: 1811-1825

[14]Mammini P,Nordt A,Holmes B,et al. Sensitivity evaluation of mounting optics using elastomer and bipod flexures [J]. Proceedings of SPIE,2003,5176: 26-35

[15]Kornmann M,Genet M. Verification of a dimensionally stable space structure [J]. Acta Astronautica,2001,48(1): 29-36

[16]Edesona R L,Aglietti G S,Tatnall A R L. Dimensional stability of materials subject to random vibration [J]. Precision Engineering,2013,37: 323-335

[17]ECSS Secretariat. ECSS-E-HB-32-26A Spacecraft mechanical loads analysis handbook [S]. Noordwijk: ESA Requirements and Standards Division, 2013: 304-325

[18]左博,范立佳,楊松,等.高分二號衛星高精度結構熱致變形分析研究[J].航天器工程,2015,24(6):64-70

Zuo Bo,Fan Lijia,Yang Song, et al. Research on high precision thermal-induced deformation of Gaofen-2 satellite structural [J]. Spacecraft Engineering,2015,24(6): 64-70 (in Chinese)

[19]盛驟,謝式千,潘承毅.概率論與數理統計[M].北京:高等教育出版社,2001

Sheng Zhou,Xie Shiqian,Pan Chengyi. Probability and statistics [M]. Beijing: Higher Education Press,2001 (in Chinese)

[20]孫振綺,張憲君.空間解析幾何與線性代數[M].北京:機械工業出版社,2011

Sun Zhenqi,Zhang Xianjun. Analytic geometry of space and linear algebra [M]. Beijing: China Machine Press,2011 (in Chinese)

(編輯:夏光)

Thermal Deformation Analysis Method of in Orbit Whole Cycle for Spacecraft

LIU Guoqing LUO Wenbo TONG Yelong FAN Lijia

(Beijing Institute of Spacecraft System Engineering, Beijing 100094,China)

Thermal deformation analysis method of in orbit whole cycle is presented. Mapping of thermal model to mechanical model during the whole cycle is executed based on method of heat conduction. Through mathematic fitting algorithm, analysis data are processed and analyzed, and process is automated by procedure. Thermal deformation analysis of whole cycle is analyzed for a remote sensing satellite. According to the analysis results, time cost of mapping is reduced from several days in the past to several hours nowadays. Besides, simulating precision is improved significantly, mapping error is less than 1%, and results can coverage cases of the whole cycle compared with the thermal deformation analysis method based on extreme load case. Dimension of thermal deformation and change rule of deformation are obtained from the calculating results meanwhile. All of the achievements are of great importance to the design of thermal stability for spacecraft.

spacecraft;in orbit whole cycle;thermal deformation;stability

2016-08-08;

2016-11-21

國家民用空間基礎設施(發改高技[2015]2429號)

劉國青,男,工程師,從事航天器機械系統設計與驗證工作。Email:liuguoqing2011@163.com。

V414

A

10.3969/j.issn.1673-8748.2016.06.007

猜你喜歡
變形結構分析
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
隱蔽失效適航要求符合性驗證分析
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
“我”的變形計
例談拼圖與整式變形
會變形的餅
電力系統及其自動化發展趨勢分析
論《日出》的結構
主站蜘蛛池模板: 亚洲人成在线精品| 欧美日韩另类在线| 1024国产在线| 无码福利视频| 国产日韩欧美视频| 国产91透明丝袜美腿在线| 992Tv视频国产精品| 波多野衣结在线精品二区| 亚洲午夜天堂| 一级毛片视频免费| 妇女自拍偷自拍亚洲精品| 日本午夜影院| 精品国产一区91在线| 欧美成人第一页| 久久夜色精品国产嚕嚕亚洲av| 精品国产成人三级在线观看| 国产成人精品一区二区| 首页亚洲国产丝袜长腿综合| 日韩午夜福利在线观看| 老司国产精品视频91| 亚洲不卡av中文在线| 亚洲熟妇AV日韩熟妇在线| 日本一本正道综合久久dvd| 亚洲精品第一页不卡| 一本大道视频精品人妻| 精品国产91爱| 色婷婷色丁香| 国产精品中文免费福利| 国产成人AV男人的天堂| 亚洲IV视频免费在线光看| 3D动漫精品啪啪一区二区下载| 欧美日韩北条麻妃一区二区| 青青草91视频| a国产精品| 日本道中文字幕久久一区| 99福利视频导航| 国产精品欧美在线观看| 久久久久88色偷偷| 国产精品无码久久久久久| 一级成人a做片免费| 中文字幕欧美日韩| 日本在线视频免费| 在线播放真实国产乱子伦| 国产在线第二页| 久久精品国产在热久久2019| 久久亚洲中文字幕精品一区| 欧美日韩在线第一页| 午夜国产大片免费观看| 国产亚洲视频免费播放| 午夜精品一区二区蜜桃| 欧美伊人色综合久久天天| 老司国产精品视频91| 中文字幕不卡免费高清视频| 亚洲无线视频| 欧美在线伊人| 女人一级毛片| …亚洲 欧洲 另类 春色| 欧美69视频在线| 波多野结衣的av一区二区三区| 欧亚日韩Av| 国产成人久视频免费| 黄色网页在线播放| 毛片在线播放网址| 精品色综合| 亚洲三级成人| 91精品啪在线观看国产| 免费又爽又刺激高潮网址| 国产亚洲精| 亚洲综合极品香蕉久久网| 91原创视频在线| 国产精品第一区在线观看| 中日韩一区二区三区中文免费视频| 精品少妇三级亚洲| 视频一区视频二区中文精品| 91小视频在线| 欧美性爱精品一区二区三区 | 久久99国产乱子伦精品免| 全色黄大色大片免费久久老太| 国产免费怡红院视频| 亚洲三级视频在线观看| 国产精品理论片| …亚洲 欧洲 另类 春色|